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A SPACE-CHARGE-FLOW COMPUTER PROGRAM 
by Carl D. Bogart and Edward A. Richley 
Lewis Research Center 


SUMMARY 

A description of a computer program for the solution of two-dimensional and axisym- 
metric space- charge-flow problems is given. The program is written in FORTRAN IV 
for an IBM 7094 computer. Primary emphasis is placed on nomenclature, word defini- 
tions, and programing procedures. The program overlay chart and subroutine flow dia- 
grams are presented, and example sets of input and output data for a variety of boundary 
conditions are given and fully explained. The program has been applied extensively for 
the analysis of various electrostatic-thrustor-design concepts, wherein positive charged- 
particle flow is conventional. It is equally applicable to various other problems, such as 
electron tube design in which the particles (electrons) have a negative charge. Output 
data listings include detailed information on potential distributions, particle trajectories, 
thrust, power and current densities, and other parameters of interest. 


INTRODUCTION 

Design requirements of electrostatic thrustors that are under investigation at the 
Lewis Research Center have given rise to the need for a numerical method of analysis of 
charged- particle trajectories subject to various boundary conditions. To fill this need, 
a computer program has been developed at Lewis that is capable of solving both two- 
dimensional and axisymmetric problems for charged-particle- flow conditions that range 
from zero charge density to space- charge limited. The program is presented herein. 

Early versions of the program have been reported previously in references 1 and 2. 
Reference 1 deals with the two-dimensional space- charge- limited- flow problem and in- 
cludes a program written in FORTRAN II for an IBM 7094 computer. In reference 2, a 
program for a special axisymmetric space-charge-limited- flow problem is given, and 
in a similar manner, other investigators have prepared computer programs for particu- 
lar space- charge- flow problems (e. g. , ref. 3). Subsequently, several additional fea- 
tures have been incorporated into the programs given in references 1 and 2, and they 



have now been merged into one overall program. 

In references 1 and 2, emphasis was placed on the mathematical techniques that were 
employed to obtain numerical solutions of the various problems and the interpretation of 
results from the physical viewpoint. Very little attention was given to a description of 
program details from the computer programing viewpoint. Because of their length and 
complexity, this has led to difficulties regarding use of these earlier programs by inde- 
pendent investigators. In this report, the purpose is to describe the present computer 
program with emphasis on preparation of input data. 


THE PROBLEM 

In this section, the mathematical model will be established from a physical model, 
and the general method of solution will be described. While much of the information in 
this section has been reported previously (refs. 1, 2, and 4), a brief description of the 
problem and general method of solution is necessary for completeness. 

The problem under consideration is either the two-dimensional or axisymmetric 
space- charge- flow problem for positively or negatively charged particles. For the in- 
vestigator of various electrostatic-thrustor designs, numerical solutions are sought that 
will provide detailed information regarding parameters, such as potential distributions 
for various boundary conditions, ion trajectories, thrust, power, and current- density 
distributions. Many of these parameters are of interest in other fields of study, such as 
electron tube design. 


Physical Model 

The problem is formulated by first defining a region of interest of a prospective 
thrustor design such as shown in figure 1. The ideal theoretical performance of this 
thrustor as well as several other thrustor configurations is given in reference 5. In the 
design shown in figure 1, ions are formed on the downstream face of the ionizer, which 
together with the focus electrode is maintained at a high positive potential with respect 
to ground. The accel electrode is at a negative potential, and the resulting field causes 
acceleration and ejection of the ions that make up the exhaust beam. Along with the 
aforementioned parameters, the minimization of ion interception on the accel electrode 
is of particular importance to a prospective thrustor design. 


Mathematical Model 

For the purposes of the numerical analysis, the region of interest, depicted in fig- 
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Ideal potential distribution 



thrustor and ideal potential distribution. 


ure 1, is laid out to scale, bounded, and overlaid with a uniform square mesh, as shown 
in figure 2. All mathematical symbols are defined in appendix A. 

Boundary conditions may be given either in terms of a potential (such as along the 
ionizer and electrodes) or as a zero normal derivative of potential. This latter specifica- 
tion is applicable, for example, along the lower boundary of symmetry and is also a rea- 
sonable approximation along the upper boundary, particularly if the boundary is normal 
to the electrode surfaces and is sufficiently removed from the space- charge- flow region. 
Thus, the best location of some of the boundaries becomes a matter of judgment. For 
example, the right boundary is not physically well defined; however, for ion-thrustor 
operation, in which neutralizers are used, it can be approximated as a straight line of 
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zero potential. The effect of the lo- 
cation of this boundary on various 
solutions has been investigated and 
found to be negligible if the distance 
from the accel electrode is about the 
same as the distance between the ion- 
izer and the accelerator (ref. 1). 

Mesh points are numbered verti- 
cally starting at the upper left corner. 
The mesh may be made as fine or as 
coarse as desired within practical 
limits. The total number of mesh 
points employed is limited only by the 
storage capacity of the computer. 
Typical problems have used from 
1500 to 3000 points. 

Additional specifications required 
for preparation of input data are the 
calculated values of the various e *s 
and 6’s (see sketch (a)), which are 
fractions of mesh- spacing distances 
measured from the mesh point under 
consideration to the boundary at a 
fixed potential. In these problem 
specifications, it is noteworthy that no restrictions are imposed with respect to the shape 
of the region. 

The equation that must be solved for the bounded region is the Poisson equation, 
which can be written for the two-dimensional problem as 



Region of 
interest 


-V 2 cp(x,y) = — p(<p,x,y) 


( 1 ) 


and for the axisymmetric problem as 

-V 2 <p(r,z) = — p(<p, r,z) 


( 2 ) 


The functions (p and p are the continuous potential distribution and space- charge- 
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density distribution functions, respectively. For the numerical, or discrete, problem, 
equations (1) and (2) are rewritten as 

-V 2 w(x, y) = f(w, x, y) (3) 


-V 2 w(r, z) = f(w, r, z) (4) 

A five-point formula approximation of equation 
(3) or (4) is made for each mesh point of the bounded 
region. For example, the equation for a simple inte- 
rior point, such as shown in sketch (b), is for the two- 
dimensional problem 

w 5 “ ^ ( w 2 + w 4 + w 6 + w g ) = ~ f 5 h 2 (5) 

For the axisymmetric problem, the equation for point 5 is 

T_ n 

r(4w 5 - w 2 - w 4 - w 6 - w g ) + - (w 6 - w 4 ) = rfgh (6) 

2 

Derivations of these equations as well as variations that arise for points near the various 
boundaries are given in references 1 and 2. 

The equations, when written for each of N mesh points, form a set of N linear al- 
gebraic equations with N unknowns. These equations, arranged in proper order, may 
be written in matrix form as 


and 
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(b) 


Aw = k (7) 

where A is a matrix consisting of the coefficients of the w’s, and w and k are column 
vectors representing the potential distribution and space- charge- density distribution 
functions. The column vector k also contains the known boundary values. The method 
used for solving equation (7) is the Cyclic Chebyshev Semi- Iterative Method, and details 
regarding the mathematical techniques used are given in references 1 and 2. 


GENERAL METHOD OF SOLUTION 

The numerical solution of equation (7) in the absence of space charge presents no 
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particular problem, and the Laplace potential distribution is readily obtained by straight- 
forward iteration. With space charge present (due to the charged particle beam), the so- 
lution is more complex since values of the space- charge function that are required in 
equation (7) are not known a priori. 

To begin the calculations that include space- charge effects, it is noted that 

f(w, x, y) = — M- (8) 

e 0 v(x,y) 

where j and v are the current-density and ion-speed functions. Similar equations ap- 
ply in axisymmetric notation. Thus, if the functions j and v can be determined, values 
of f may be calculated for use in equation (7). The j and v functions are not known 
explicitly; however, they may be initially approximated from known physical laws, and 
additional iteration techniques may be employed. The initial j and v estimates are ob- 
tained with the aid of the Laplace potential distribution. 

As discussed in reference 4, solutions may be obtained for either space-charge- 
limited-flow or less-than-space-charge-limited-flow problems. 


Space -Charge -Limited Problems 

The procedure to obtain values of j and v is as follows: The ionizer is divided into 
an arbitrary number of line segments of length s, as shown in sketch (c). Ion trajecto- 
ries, which form the boundaries of 
stream tubes, are assumed to start 
from the ends of each line segment. 

As the trajectory moves to the right to 
each mesh column, its position is first 
estimated from the slope given by the 
velocity components at the previous 
z mesh column. Its final position is de- 
termined by an iterative procedure that 
employs the equations of motion. The procedure and various trajectory possibilities are 
described in detail in reference 1. Velocities at any point in the region are determined 
from the law of conservation of energy. 

The current in each stream tube is found at the ionizer (see sketch (c)) by a parallel- 
plate approximation where the "plates" are assumed separated by a distance £ and a 
potential difference of A cp. From Child’s law, the space- charge-limited current density 
is given as 
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(9) 


j 


E “ 


4 l2 q A <p_ 
9 °\m £ 2 


In this calculation, either A <p or S. may be held constant for the various stream 
tubes. And as in the trajectory calculation, the Laplace potential distribution is used to 

■ obtain the initial A cp and £ values. 
Equation (9) , if multiplied by the 
stream-tube cross-sectional area, 
gives the total current flowing in each 
stream tube. With the stream-tube 
current determined, the j (x, y) distri- 
bution is calculated as follows: Con- 
sider a typical mesh point in the region 
* of interest as shown in sketch (d). To 
find the space- charge function fg 
chargeable to the subregion associated with point 5, it is necessary to sum the current 
contributions from each tube and divide by the subregion cross-sectional area (h x unit 
depth) and the magnitude of the average ion velocity in the x-direction in each stream 
tube. Thus, referring to sketch (d) gives 



Subregion of 
influence 
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J a'b ) J bc 

ft-).* ( v 4c 



( 10 ) 


where 



and where (v ) , (v„) , and (v_ ) are the x-components of the average ion velocity 

\ x/ a , b \ x/ bc \ x/ cd . 

and ab, a’b, be, etc., are the stream-tube and/or line-segment designations, as shown 
in sketch (d). 
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The current density of the subregion is calculated by dividing by the subregion cross- 
sectional area taken normal to the x-direction; therefore, the x-components of the speeds 
must be used in equation (10). 

It is also convenient in the axisymmetric problem to approximate the current density 
at the ionizer by using equation (9). The total current in each stream tube, however, 

must now be calculated as J t = jg 27 rrs, 
where r is the average radius to the ionizer 
line segment. For the space- charge- density 
calculation, area ratios must be used rather 
than the simple line- segment ratios employed 
in the two-dimensional problem. For exam- 
ple, consider the typical mesh point shown in 
sketch (e). Here, 



where 



( 12 ) 


(13) 


so that 



(14) 


To recapitulate, space- charge- limited solutions to equation (7) are obtained by an 
iteration procedure in which the Laplace potential distribution is used to calculate initial 
ion trajectories and to obtain initial estimates of the current flowing in the stream tubes. 
The stream-tube current density at the ionizer is calculated from a Child’s law relation 
for space-charge-limited-current flow. Values of the space- charge-density distribution 
are then determined and supplied in equation (7), which is solved to provide a new poten- 
tial distribution. This completes one cycle. The new potential distribution then forms a 
basis for recalculation of trajectories and current densities, and the process continues 
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until convergence is obtained. Details of the factors that affect the rate of convergence 
are given in reference 1. In general, from five to seven cycles are required with a total 
machine time of about 5 minutes. 


Problems Less Than Space-Charge Limited 

Previously, the current-density distribution of equation (8) was determined by first 
applying equation (9) to each stream tube and then using equation (10) or (14). Each sub- 
sequent iteration cycle of equation (7) produced a new potential distribution that is used 
to calculate new values of the emitter current density jg. 

On the other hand, jg may be prescribed at some value, less than the space-charge- 
limited value, and held constant throughout the iterative procedure. This implies that 
the current in each stream tube remains constant. The current-density-distribution and 
velocity- distribution calculations required in equation (8) may then proceed in the same 
manner as described. 


COMPUTER PROGRAM 


LINKO 

I 

MATRIX 


r r 


rfir 


The main program, LINKO, is divided into three major core loads, as shown in fig- 
ure 3; each core load consists of several subroutines. Core load 1 is for data input and 
initialization. Core load 2 contains the subroutines used for calculation of the potential 
distribution. The charged-particle trajectory coordinates and space-charge-density dis- 
tributions are calculated in core 
load 3. 

Common statement symbols 
used in the program are defined in 
appendix B, and a complete listing 
of the program is given in appen- 
dix C. Brief descriptions of each 
of the core loads and subroutines 
along with schematic flow diagrams 
are presented in appendix D. An 
understanding of the input data prepa- 
ration and output data listings is es- 
sential to use the program, and 
these topics will be discussed in de- 
tail in the sections that follow. 
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Figure 3. - Overlay chart of program. 
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Figure 4. - Model for axisymmetric sample problem. Mesh spacing, 0.25 units; focuser and emitter potential, 1 kilovolt; accelerator potential, 
-1.0 kilovolt; right boundary, 0 kilovolt; normal derivative, zero along skew boundary and upper and lower boundaries. 


Input Data 

In this section the quantities required for preparation of input data are described. 
The explanation is based on the relatively simple model, shown in figure 4, that will be 
used as a sample problem for results to be presented in the next section. 

In the explanation that follows, various subroutines appearing in core load 1 (fig. 3) 
that require input data are presented. Words are presented in the order in which they 
appear in each subroutine. Mesh point and mesh column numbering should always be or- 
dered as shown in figure 4. The card column locations for the sample input data will be 
indicated in the input data listing. 

It is important to note that, although the x, y or r, z coordinate orientation is as 
shown in figure 4, the program works only with x, y coordinates. Thus, unless other- 
wise noted in the following discussion, data input is always in terms of x, y coordinates. 


Subroutine CHN14 

NPIT. - If a first guess of the potential distribution is to be read in from subroutine 
GUESS, set NPIT = 0. If the potential distribution is available from subroutine BCDUMP 
(from prior computer run), set NPIT = 1. 

NTP. - The absolute value of the maximum JT-type number plus 4 is given by NTP; 
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that is, NTP = (j JT max I + 4). The JT-type numbers are defined in the subsection sub- 
routine XTCAL. 

KAT - - This is a mesh column number and is used in connection with the test for cur- 
rent impingement on the first electrode for which an impingement test is desired (e. g. , 
the accel electrode in fig. 4, p. 10). Of the mesh columns passing through the electrode, 
KAT is the minimum mesh column number where impingement can occur. If no impinge- 
ment test is desired, set KAT = 0. 

KATT . - Similar to KAT, KATT is the maximum mesh column number in the first 
electrode for which an impingement test is desired. If no impingement test is desired, 
set KATT = 0. 

NT . - The total number of mesh points is given by NT. 

NTJ . - Trajectories are equally spaced along the emitter (a spacing equivalent to 
1 mesh width is usually adequate), and NTJ is equal to the number of spaces. 

NT A . - Number of pairs of x, y coordinates used to specify the emitter surface is 
given by NTA. For the model shown in figure 4, only the beginning and the end pairs are 
needed since the emitter is a straight line. 

KAN . - KAN is a mesh point number ^ to 1^ mesh widths away from the emitter sur- 
face centrally located with respect to the emitter surface. The space-charge-density 
function is checked for convergence at this point. 

KB A . - KBA is a mesh point number such that if an equipotential line is taken through 
this mesh point, the equipotential line will be no closer than 1 mesh width from the emit- 
ter and no farther than 3 mesh widths. It may not be possible to satisfy these conditions 
simultaneously for some configurations in which case the selection of KBA should be 
based on the former condition. This equipotential line is used to calculate the emitter 
current density. 

KAB . - KAB is a mesh column number. It is chosen to be the farthest mesh column 
to the right of the emitter required to establish a region in which the equipotential line 
will be located that is used to calculate the emitter current density. 

KAT1 . - Similar to KAT, KAT1 is the smallest mesh column number in the second 
electrode for the impingement test (if there is only one electrode as in fig. 4, or, if no 
impingement test is desired on the second electrode, set KAT 1 = 0). 

KAT2 . - Similar to KATT, KAT2 is the largest mesh column number in the second 
electrode for the impingement test (if there is only one electrode, or if no impingement 
test is desired on the second electrode, set KAT2 = 0). 

IAS. - For two-dimensional problems, set IAS = 0. For axisymmetric problems, 
set IAS = 1. 

NPL. - The number of mesh points per mesh column is given by NPL. 

KBB . - To calculate space-charge-limited current density, set KBB = 0. To specify 
the current density along the emitter set KBB = 1. At this point, five heading cards 
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(card columns 1 to 72) are read into the program and printed out at the beginning of the 
data output listing. 

JAS(J) . - The Laplace potential distribution leads to an over estimation of the space- 
charge density and subsequently to a potential distribution that may be over -space -charge 
limited (i. e. , there are potentials at mesh points near the emitter that are higher than 
the emitter potential) . To check this condition, it is necessary only to examine a few 
points near the emitter. Thus, JAS(l) is the number of mesh points to be checked and 
JAS(2) to JAS[JAS(1) +1] are the mesh point numbers. A maximum of 19 points can be 
checked with the present dimension of JAS. 

AW . - Atomic weight (amu) of charged particles under consideration is given by AW, 
which is positive for positively charged particles and negative for negatively charged 
particles. 

VA . - The emitter potential in volts is given by VA. 

H. - The mesh size H is a fraction of chosen unit of length. Results will be in terms 
of whatever dimension is used for H. 

VAT . - The potential of the maximum equipotential line to be printed out in volts is 
given by VAT. 

VBT . - The potential of the minimum equipotential line to be printed out in volts is 
given by VBT. 

SIZE . - This is the step size of equipotential lines to be printed out in volts. If 
SIZE = 0, no equipotentials are printed out. 

RCU(J) . - Read in only when the current density is specified (KBB = 1) along the 
emitter surface, that is, less-than-space-charge- limited problems. The current density 
at the emitter of the J in segment of the emitter is RCU(J). The segments are bounded by 
charged- particle trajectories and are numbered consecutively, starting from 1, and in- 
crease with increasing y (decreasing r). Note for example, that if H is given in milli- 
meters, the RCU(J)'s must be specified as amperes per millimeter squared. 

At this point, the data from subroutine XTCAL is read in. 

ATX(J) . - This is the x-coordinate of the pair of x, y coordinates given to 
specify the emitter surface. 

ATY(J) . - This is the y-coordinate of the J Ln pair of x, y coordinates given to 
specify the emitter surface. By convention, ATY(J) < ATY(J + 1) and usually these are 
coordinates of y-mesh lines except possibly for the first and/or last values. It should be 
emphasized that for axisymmetric problems it is the y-coordinate that is used (not the 
r- coordinate). 

ER(J) . - A set of y- coordinates of the electrode edges are required to calculate im- 
pingement on the electrodes. The largest y-coordinate of the first mesh column in the 
first electrode is ER(1) (corresponds to mesh column KAT) . Similarly, ER(2) to 
ER(KATT - KAT + 1) are the largest y-coordinates for the remaining mesh columns in 
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the first electrode. Likewise ER(KATT - KAT + 2) to ER(KAT2 + KATT + 2 - KAT - KATl) 
are the y- coordinates for impingement on the second electrode, should there be a second 
electrode. 

The data from subroutine JTINT is read in after the ER's. 

XR . - This is the spectral radius of iteration matrix and it is the absolute value of 
the largest eigenvalue of the iteration matrix (see ref. 1). It is generated in subroutine 
EVC and is printed as part of the data output. Its value depends only on the physical con- 
figuration. Initially, a blank card is included to satisfy the read statement. For subse- 
quent runs of the same configuration, XR is available from the data output. 


Subroutine XTCAL 

The finite-difference equation replacing Poisson's equation at the mesh point N has 
the form 

U(N) = C 1 (N)U(N + Jj) + C 2 (N)U(N + J 2 ) + C 3 (N)U(N + Jg) + C 4 (N)U(N + J 4 ) + C 5 (N)RH(N) 

(15) 

where the U(N + J.), i = 1, • • • , 4, are potential values, the C^(N), i = 1, • • • , 5, are 
coefficient weights, and RH(N) is the charge density at the N n point. Subroutine XTCAL 
calculates the coefficients in equation (15), C j (N) , • • •, Cg(N), for the different mesh- 
point configurations. Generally, the coefficients depend on the distance from the mesh 
point to a boundary at a fixed potential and, in the axisymmetric case, on the radius as 
well. For either the two-dimensional or axisymmetric problem, the data input has the 
following format for each card: 

KT(1), KT(2), KT(3), KT(4), XF1, XF2, XF3, XF4, R. 

KT(1) to KT(4) are relative numbers that are added to the mesh point number. They are 
the numbers J t , • • ■ , J 4 in equation (15). XF1 to XF4 are the positive distances 

to either the boundary at a fixed potential or the nearest 
mesh point. The distance to the central mesh point from 
the axis of symmetry in the axisymmetric case is R (no 
entry required in two-dimensional problems). These 
quantities are further illustrated in sketch (f). 

More specifically, KT(1) is used to select the upper 
vertical mesh point, and its value is usually equal to -1. 
For example, from examining figure 4 (p. 10), it is seen 
that for all mesh-point configurations, except for the ones 


X ** 

from the central mesh point N 
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along the upper boundary and the one associated with mesh point 24, KT(1) = -1. For 
mesh point 24, it is necessary to use an entry for KT(1). that will select a mesh point at a 
fixed potential (i. e. , a mesh point located inside the accel electrode). Thus, KT(1) can 
have the entry 3 which picks out point 27, or 7, which picks out point 31. The XFl's are 
related to the KT(l)’s, and since they are zero for the upper boundary mesh configura- 
tions, the values used for these KT's are not important. It is convenient to use -1. A 
similar argument applies to other boundaries. The relative number used to select the 
lower vertical mesh point to be used in equation (15) is KT(2), and examination of fig- 
ure 4 reveals that KT(2) = 1 except for mesh point 23. At mesh point 23, KT(2) can either 
be 4 or 8. The relative number picking out the horizontal point to the right is KT(3), and, 
in general, it is merely equal to the number of points per column. In every case in fig- 
ure 4, KT(3) is 5 except for mesh point 23. For mesh point 23, KT(3) can be either 4 
or 8, since as before, it is only necessary to pick out some point in the accel electrode 
so that the correct potential value will be used in equation (15). The relative number 
picking out the horizontal point to the left is KT(4), and, in general, it is equal to minus 
the number of points per column. In figure 4, KT(4) is seen to be -5 in every instance 
except for point 28 in which case it can either be -1 or 3. It should be mentioned that 
while the preceding remarks pertaining to mesh points 23 and 24 may make it appear as 
if these two points are rather special, an important feature illustrating the flexibility of 
the program has been demonstrated; that is, with this method of input data preparation, 
electrodes may be "infinitely thin, " the only requirement being that a mesh point at the 
desired potential must exist somewhere in the array. 

The distances XFl, • • •, XF4, are measured, respectively, in the upward vertical 
(-y) direction, the right horizontal (+x) direction, the downward vertical (+y) direction, 
and, the left horizontal (-x) direction (see sketch (f)). If any of the XF’s equal zero, it 
is assumed in the program that the normal derivative equal to zero is specified in the 
corresponding direction, and the appropriate formula (refs. 1 and 2) is then used to cal- 
culate the coefficients. When it is desired to have a skew boundary that has the normal 
derivative equal to zero as the "boundary condition, " the pertinent boundary points 
are described somewhat differently. Sketch (g) depicts a portion of the region from fig- 
ure 4 that contains a skew boundary. Note in 
sketch (g) that a normal to line AB from 
point 11 intersects between points 7 and 12 
at point X, which is at a distance DA from 
point 12 and a distance DB from point 7. 
Similarly, a normal to AB drawn from 
B point 16 intersects at Y, which is at a dis- 

tance DA from point 17 and DB from point 12. 
The normal derivative equal to zero 
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(dU/dn = 0) implies that an equipotential line passing through points 11 or 16 would be at 
right angles to the line AB and that if the potential at point 11, for example, were extrapo- 
lated along the line 11 X, the potential at point X would be at the same value. The poten- 
tial at point 11 is not known but the potential at point X can be found by interpolation be- 
tween values at points 7 and 12. Values of potential at points 7 and 12 are available be- 
cause they are "interior points. " Similarly, the potential at point 16 is the same as that 
at point Y. Thus, using linear interpolation gives the potential at point 11 as 

U(U) = (—^ V(7) + f—® 1 * Vj(12) (17) 

\DA + DB/ \DA + DB / 

and at point 16 as 

U(16) = ( DA V( 12 ) +( — — )U(T7) (18) 

\DA + DB/ VDA + DB/ 

Recall that the general form of equation (15) is U(N) = C^(N)U(N + Jj) + C 2 (N)U(N+ J 2 ) 
+ • • • + C 5 (N)RH(N). The C^N), in general, are calculated from the finite-difference 
equations (by using the XF distances), but it is a desirable feature to be able to enter 
certain C^(N) coefficients (e. g. , those associated with eqs. (17) and (18)) directly into 
the computer. This feature is accomplished by setting XF1 = -C^(N), XF2 = CgW, 

XF3 = C 3 (N), XF4 = C 4 (N), and R = C 5 (N). The XT's of this card can then be assigned 
relative numbers to associate the coefficients with the corresponding mesh points, that 
is, for these special cards, KT(1) is associated with XF1, KT(2) with XF2, etc. Refer- 
ring to equations (17) and (18) shows that the entries for the mesh configuration of mesh 

points 11 and 16 would then have the form -4, 1, 0, 0, -( — V ( — 5® V 0, 0, 0. 

\PA + DB / \DA + DB/ 

Note that mesh points 12 and 17 are considered as regular points since the distances XF1 
to XF4 are only different from a mesh spacing when measured to a boundary where the 
potential is held constant or in the special cases just described. 

It will become evident in proceeding from mesh point to mesh point that duplication 
of all nine entries on each card can occur. To eliminate this duplication, all similar 
mesh points are assigned a JT-type number of the form ±(5n + 4), n = 1, 2, • • • , (see 
fig. 4, p. 10). In this manner, a large number of mesh points can be represented by one 
data card. The method used to get the JT-type numbers into the program as data is de- 
scribed in the subroutine JTINT section. 

As is shown in figure 4, the assignment of actual JT-type numbers is arbitrary as 
far as correlation with mesh point numbers. However, in correlating XTCAL data cards, 
the approach is to associate the first data card for XTCAL with mesh points of type ±9, 
the next, type ±14, etc. No XTCAL data cards are required for mesh points held at a 
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fixed potential, although these points are also assigned JT-type numbers. The numbers 
0 to 6 are reserved for these points. All points associated with the focus electrode may 
be assigned the JT-type number 0; all accel electrode points, type number 1; and so on 
to a maximum of type number 6. These mesh points will then have the potential values 
as prescribed in subroutine GUESS. 

The JT-type numbers are of the form ±(5n + 4) because there are exactly five coeffi- 
cients (eq. (1 5)) necessary for each mesh-point configuration. The sign designation, 
plus or minus, is used to describe the central mesh points further. A positive type num- 
ber is used to designate mesh points whose coefficients C^(N + J.) are calculated from 
the finite-difference equations. A negative type number is used to indicate those points 
that are free to change in value but have coefficients calculated from linear interpolation. 
Examples of this type of point are shown in figure 4 (p. 10), that is, points 11, 16, 21, 

22, and 26. As previously discussed, they arise because of the skew boundary that has 
the boundary condition of normal derivative equal to zero specified. It is not necessary 
to make consecutive duplicate data cards: only one card per different mesh- point configu- 
ration is required. Also, no card is required for mesh points held at a fixed potential. 


Subroutine JTINT 

Once a JT-type number is assigned to each mesh point, the computer generates the 
matrix equation 


IU = CU + K (16) 

where, if M denotes the number of mesh points, I is the M x M identity matrix, C is 
the M x M matrix having entries in each row Cj(N), C 2 (N), Cg(N), and C^(N); K is the 
1 x M column vector with entries Cg(N)RH(N), and U is the lxM column vector of 
potentials. If the C matrix contained M x M nonzero entries, it would be prohibitive to 
store C in core, but C contains only four entries in each row (see eq. (15)) so that the 
C matrix can be specified by at most 4M numbers. 

To enter the JT-type numbers into the program as data, they must be punched on 
cards. While it is possible to write out the JT-type numbers, one for each mesh point, 
punch them on cards, and read them into the JT-array, this method is inefficient for 
large M. A more general method is used. 

Description o f the JTINT dat a card s. - The first card has only one entry, JA, which 
is the number of cards that contain data for the JT array. The JT cards all have the 
same format, that is, KA; KB(1) to KB(13). For the axisymmetric geometry, KA is 
always equal to 1; the two-dimensional variation of KA will be discussed later. The 
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number of consecutive mesh points to be assigned the JT-type number in KB(2g) is 
KB(2g - 1), g= I,- • 6. 

The convention for setting up the KB array for an axisymmetric problem is to start 
at mesh point number 1 and count the number of consecutive mesh points in the first row 
(not column) that have the same JT-type number. Looking at the first row in figure 4 
shows that this number, KB(1), is 2 since points 1 and 6 are boundary points at a fixed 
potential. Thus, the entry of KB(1) is 2, and KB (2) has the entry that is the JT-type num- 
ber assigned to mesh points that are held at the focuser potential. Then KB(3) is seen to 
be 4, and KB(4) is the JT-type number assigned to these mesh points; KB(5) is 2, and 
KB(6) is the JT-type number assigned to these mesh points; KB(7) is 3, and KB(8) is the 
JT-type number associated with points 41, 46, and 51; KB(9) is 1, and KB(10) is the 
JT-type number associated with point 56. This completes the specification of the points 
of row 1; however, note that KB(ll) and KB(12) have not yet been assigned. The specifi- 
cations for row 2 may begin in KB(ll) and KB(12), or they may be left blank and row 2 
specifications begun on the next card with KB(1) and KB(2). The procedure is continued 
until all mesh points have been covered. Recall that the number of cards generated in 
this process is JA. 

For an axisymmetric problem (IAS =1), it is also necessary to supply a final card 
giving the distance from the axis of symmetry of the model to the top of the region. This 
number is denoted in the program by BASE. 

For two-dimensional problems, the process is similar except that it starts at mesh 
point 1 and goes down the first column in assigning values to KB(K), K = 1, • • • , 12, 
then down the second column, etc. The reason for going down the columns in two- 
dimensional problems is that it is frequently possible to repeat the KB specifications of 
a card several times. 

Rather than punching the same card over and over again, the number of repetitions 
desired may be specified by KA. Note that utilization of this feature is not possible for 
axisymmetric problems wherein the coefficients of equation (15) vary with respect to 
radius. Thus, if a columnwise procedure were used, each mesh point would require in- 
dividual specification. JT-type number assignments for a two-dimensional problem are 
shown in figure 5. 


Subroutine GUESS 

Subroutine GUESS is used to "initialize’ * the potential field. Data are read in only 
if NPIT = 0. 

GEP(J) . - As discussed in the data input section for subroutine XTCAL, boundary 
points that are held at a specified potential are assigned JT-type numbers from 0 to 6. 
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Mesh points designated as type 0 are assigned the value of GEP(l), those of type 1 the 
value of GEP(2), etc. All other mesh points are initially given the value zero. 


SAMPLE PROBLEMS 

In the previous section, the input data words were presented and described. In this 
section a sample problem based on the axisymmetric model shown in figure 4 will be de- 
scribed; that is, numerical values will be assigned to the input data words, and results 
will be given. Machine listings of the input and output are presented in appendix E. 

The mesh size shown in figure 4 is far too coarse to be used for meaningful analysis 
in the physical sense; however, for purposes of illustration of data preparation it is ap- 
propriate. As a result of the coarse mesh, the accelerator tip (shown as the shaded re- 
gion in fig. 4, p. 10) extends into a mesh square and is not accounted for in impingement 
calculations. It will be assumed first that the model represents an electrostatic thrustor, 
which uses cesium as a propellant, and that the current density will be specified along 
the emitter. Voltages and distances will be as shown in figure 4. A set of input and out- 
put data for the same model treated as a space- charge- limited two-dimensional problem 
is also included in appendix E. 


Input Data Preparation 

The numerical values of the input quantities are given next. The data are presented 
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in the order in which they are read in by the program (i. e. , the order in which they ap- 


pear in appendix E) . 


Subroutine CHN14. - 


NPIT 

0 

A first guess of the potential distribution will be read in from sub- 
routine GUESS. 

NTP 

78 

(| JT max |+4). The numerical value of NTP cannot be assigned 
until after the JT-type numbers have been determined (see 
fig. 4, p. 10). 

KAT 

4 

The first column for impingement test 

KATT 

6 

The last column through the accel electrode for which an impinge- 
ment test is desired 

NT 

60 

Total number of mesh points 

NTJ 

4 

Number of spaces between trajectories 

NTA 

2 

Because the emitter is a straight vertical line, only two pairs of 
coordinates are required to specify it. 

KAN 

9 

Mesh point for convergence check 

KBA 

10 

Mesh point locating equipotential for current density calculation 

KAB 

2 

Mesh column locating region for KBA equipotential 

KAT1 

0 

No second electrode for impingement 

KAT 2 

0 

No second electrode for impingement 

IAS 

1 

Designates axisymmetric problem 

NPL 

5 

Number of mesh points per column 

KBB 

1 

Current density will be specified. 

At this point, 

five heading cards are read. 

JAS(l) 

3 

Number of mesh points to check for over space charge 

JAS(2) 

8 

Mesh point to check for over space charge 

JAS(3) 

9 

Mesh point to check for over space charge 

JAS(4) 

JAS(5) 

10 

Mesh point to check for over space charge 

to > 

JAS(20) J 

0 

Only three mesh points are checked for over space charge limited 
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Refer to appendix E and note that two cards are required to satisfy the READ statement 
for the JAS-array: 


AW 

132. 91 

amu 

VA 

1000. 

V 

H 

0.25 

arbitrary units 

VAT 

1000. 

V 

VBT 

-1000. 

V 

SIZE 

200. 

V 

RCU(l) 

1. 0 E-4 

A/unit^ 

RCU(2) 

1.0 E-4 

A/unit ^ 

RCU(3) 

1. 0 E-4 

A/unit ^ 

RCU(4) 

1.0 E-4 

A/unit^ 


As previously discussed in assigning JT-type numbers, the manner in which they are 
associated with the mesh points is arbitrary, the only requirement being that mesh points 

which are similar must be assigned the same type number. In 
table I the JT-type numbers and associated mesh points are 
indicated. The JT-type numbers and associated mesh points 
for which data are not necessary are given in the table at the 
left. For reference, the JT-type numbers may be punched on 

TABLE I. - XTCAL DATA CARDS 


JT-type number 

Mesh points 

0 

1 to 6 

1 

27,31,36 

2 

56 to 60 


Data 

JT-type 

Reference mesh points 
(from fig. 4, p. 10) 

KT(1) 

(2) 

(3) 

(4) 

XFl 

XF2 

XF3 

XF4 

R 

number 

-1 

1 

5 

-5 

0.00 

.25 

.25 

.25 

1.00 

9 

41,46, 51 

-1 

1 

5 

-5 

.25 

.25 

.25 

.25 

.75 

14 

12,17,37,42,47, 52 

-1 

1 

5 

-5 

.25 

.25 

.25 

.25 

. 5 

19 

13,18,33,38,43,48, 53 

-1 

1 

5 

-5 

.25 

.25 

.25 

.25 

.25 

24 

14, 19,29, 34,39,44,49, 54 

-1 

1 

5 

-5 

.25 

.25 

0.00 

.25 

0.00 

29 

15,20,25,30,35, 40,45, 50,55 

-1 

1 

5 

-5 

.25 

.25 

.25 

.20 

. 5 

34 

8 

-1 

1 

5 

-5 

.25 

.25 

.25 

.20 

.25 

39 

9 

-1 

1 

5 

-5 

.25 

.25 

0.00 

.20 

0.00 

44 

10 

-1 

4 

4 

-5 

|| 

.25 

.025 

.050 

.25 

.5 

49 

23 

3 

1 

5 

-5 

. 125 

.25 

.25 

.25 

.25 

54 

24 

-1 

1 

5 

-1 

.0834 

.25 

.25 

.10 

. 5 

59 

28 

-1 

1 

5 

-5 

.042 

.25 

.25 

. 05 

.75 

64 

32 

-4 

1 

5 

-5 

-. 50 

. 50 

0.00 

0.00 

0.00 

-69 

11,16,21,22,26 

-1 

1 

5 

-5 

. 165 

.25 

.25 

.08 

.75 

74 

7 
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the XTCAL data cards in card columns 73 to 80. 

Subroutine CHN14 data (cont.) . - 

ATX(l), ATX(2) . 05, . 05 x- emitter coordinates 

ATY(l), ATY(2) .50,1.0 y- emitter coordinates 

ER(1), • . . , ER(3) . 625, . 400, . 23 impingement coordinates 

The data for subroutine JTINT are to be included at this point. Since the sample 
problem is being considered as an axisymmetric configuration, it will be recalled from 
the previous section that the ordering of the KB numbers will be from left to right, that 
is, along rows. 

JA 7 This value cannot be assigned until the 

following JT cards (table H) have been 
established. 

TABLE H. - JT DATA CARDS 


KA 

KB 


a) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 

(ii) 

(12) 

1 

2 

0 

4 

-69 

2 

1 

3 

9 

1 

2 



1 

1 

0 

1 

74 

2 

14 

1 

-69 

1 

1 

1 

64 

1 

4 

14 

1 

2 









1 

1 

0 

1 

34 

2 

19 

1 

49 

1 

59 

5 

19 

1 

1 

2 











1 

1 

0 

1 

39 

2 

24 

1 

54 

6 

24 

1 

2 

1 

1 

0 

1 

44 

9 

29 

1 

2 






BASE 1. 0 BASE is read in only if IAS = 1, that is, 

for axisymmetric problems. 

Subroutine CHN14 data (cont.) . - 

XR No value for the spectral radius is available, 

but a blank card is necessary here. 

Subroutine GUESS data. - 


GEP(l) 

1000. 

V 

GEP(2) 

-1000. 

V 

GEP(3) 

0 . 

V 


The preceding values represent the potential values assigned, respectively, to the 
emitter and focuser, the accel electrode, and the straight boundary at the far right of the 
region. The values are associated with JT-type numbers 0, 1, and 2. After the Laplace 
potential distribution is determined, it is punched on cards as part of the output of the 
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Figure 6. - Axisymmetric model. Focuser and emitter potential, 1 kilovolt; accelerator potential, -1 kilovolt; emitter cur- 
rent density specified at 12.7 milliamperes per square centimeter; percent impingement current at accelerator, 27. 


program. If it is desired to run the problem again with the Laplace potential as the ini- 
tial distribution, these cards should be included at this point, rather than data for sub- 
routine GUESS, and NPIT set equal to 1. 


Output Data Interpretation 

A listing of the output data for the sample problem is given in appendix E, and re- 
sults are plotted in figure 6. The listing begins with the five heading cards. The next 
portion of the data consists of the KT and XT from subroutine XTCAL. The XT 
are the coefficients of equation (15) that were calculated from the XF. Following the 
KT, XT printout is the number of iterations required to converge on the spectral radius 
XR and its value. The Laplace potential field convergence information and mesh point 
numbers along with their corresponding potential values are listed next, followed by the 
listing of the x, y coordinates of various Laplace equipotentials. 

The quantities given next are self-explanatory. It should be noted that the trajecto- 
ries listed at this point are calculated from the Laplace potential distribution. Conver- 
gence information is then given relative to the Poisson solution, where RHLOW, RHUP, 
and RH refer to the space- charge-density function at a predetermined "test" mesh point, 

KAN. The RH values at the various mesh points are the space- charge- density values 
o 

multiplied by H /4, where H is the mesh spacing. Cycles 1, 2, 3, • • * , n refer to 
successive iteration cycles, as described in reference 1. Finally, the converged Poisson 
solution is listed. 
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CONCLUDING REMARKS 


The purpose of this report has been to describe in detail a computer program capable 
of solving a wide variety of space- charge- flow problems. The approach taken has been 
one in which emphasis has not been on the mathematical relations or physical interpreta- 
tion, but rather to cover all aspects of input data preparation and output data interpreta- 
tion. Toward this end, input words were fully defined and flow diagrams presented. 
Sample problems were then used to explain the program further. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 18, 1965. 
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APPENDIX A 


MATHEMATICAL SYMBOLS 

[The units used are the International System or SI.] 


A 

matrix of eq. (7) 

w 

potential distribution function for 

f 

space-charge-density distribution 


discrete case 


function for discrete case 

w 

column vector of eq. (7) 

h 

mesh spacing 

x,y,z 

variables 

J 

current 

5,e 

fraction of mesh spacing (see 

j 

current density 


sketch (a)) 

k 

column vector of eq. (7) 

e o 

permittivity of free space 

a 

distance, defined in eq. (9) and 

P 

space- charge- density distribution 


sketch (c) 


function for continuous case 

m 

particle mass 

<P 

potential distribution function for 



continuous case 

q 

unit charge 





Subscripts: 

r 

variable 

A 

accelerator 

r 

average radius to ionizer line 




segment 

E 

ionizer 

s 

length of line segment 

t 

stream tube 

V 

ion speed 

x, z 

direction 
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APPENDIX B 


AREM 

ATX(J) 

ATY(J) 

AX 

AY(J) 

BASE 

CU(J) 

CUD(J) 

DC 

DCC 

DELY 

DX 

EPS 

ER(J) 

ETX(J) 

ETY(J) 

H 

IAS 

JAS(J) 

JD 

JOT 

JT(J) 

KAB 

KABB 

KAN 

KAT 


COMMON STATEMENT SYMBOLS 

arc distance between trajectories at emitter surface 

problem specification (see data input discussion) 

problem specification (see data input discussion) 

x- coordinate of trajectories as they are calculated 

y- coordinates of trajectories 

see data input discussion 

currents for stream tubes 

current densities at emitter surface 

distance used to sum tube currents 

distance used to sum tube currents 

mesh size in y-direction 

mesh size in x-direction 

convergence criterion for potential field 

see data input discussion 

beginning x-coordinates of trajectories at emitter surface 
beginning y-coordinates of trajectories at emitter surface 
mesh size 

see data input discussion 
see data input discussion 
intermediate storage 

printout counter for trajectory coordinates 
vector of type numbers 
see data input discussion 
see data input discussion 

problem specification (see data input discussion) 
program control word (see data input discussion) 
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KATT 

KAT1 

KAT2 

KB(J) 

KBA 

KBB 

KBF 

KCH(J) 

KRL 

KT(J) 

LAST 

LB(J) 

LC(J) 

MO 

NAJ 

NPIT 

NRD 

NRL 

NT 

NTJ 

NURL 

NXEP 

PTX(J) 

PTY(J) 

RCU(J) 

RH(J) 


program control word (see data input discussion) 
see data input discussion 
see data input discussion 

used in subroutine JTINT as problem specification (see data input discussion) 
and then as cycle print control 

see data input discussion 

see data input discussion 

indicates whether or not emitter surface extends to top boundary of region 
trajectory reflection counter 
cycle counter 

problem specification (see data input discussion) 

indicates whether or not emitter surface extends to lower boundary of region 

internal control parameters 

internal control parameters 

indicates upper or lower bound test for RH 

number of tubes 

program control word (see data input discussion) 
printout counter for trajectory coordinates 
maximum number of cycles to converge to Poisson solution 
total number of mesh points 
number of trajectories 

number of iterations on potential distribution 

mesh point number where maximum potential change is occurring 

x- coordinates of equipotential line used to calculate current density at emitter 
surface 

y- coordinates of equipotential line used to calculate current density at emitter 
surface 

see data input discussion 

present space- charge-density function is stored in RH array except in sub- 
routine EVC where an intermediate iterate on spectral radius is stored 
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RHDOWN 


RHUP 

RIN 

RX 

SAU(J) 

SEM 

SIZE 

U(J) 

UB(J) 

URH(J) 

VA 

VAT 

VBT 

VX(J) 

VY(J) 

XD 

XEP 

XK 

XMP(J) 

XQM 

XR 

XT(J) 

YEP 


lower bound on RH 
upper bound on RH 
width of region 
suppression factor for RH 
current at emitter surface 
total emitter length 

problem specification (see data input discussion) 

present potential field is stored in U array except in subroutine EVC where 
an intermediate iterate on spectral radius is stored 

potential field one cycle back is stored in UB array except in subroutine 
EVC where an intermediate iterate on spectral radius is stored 

space-charge-density function one cycle back is stored in URH array 

emitter potential 

problem specification (see data input discussion) 
problem specification (see data input discussion) 

x-velocity components of trajectories except in subroutine EQLINE where 
x- coordinates of equipotential lines are stored 

y- velocity components of trajectories except in subroutine EQLINE where 
y-coordinates of equipotential lines are stored 

intermediate storage 

storage for maximum potential change from iteration to iteration 
intermediate storage 
error factor 
charge-to-mass ratio 

problem specification (see data input discussion) 
coefficients in eq. (15) 
permittivity of free space 
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APPENDIX C 


FORTRAN LISTING 


C LINKO IS USED AS THE MAIN PROGRAM 

COMMON AREM, ATXI 51 ) ,A TY < 51 > , AX, A Y C 51 ) , BASE, CU(51), CUD(51), 

1 DC, OCC, DELY, DX, EPS, ER<40), E TX ( 5 1 ) , ET Y ( 51) , H , I AS , J AS ( 20 > , J D 

2 , JOT , J T I 3000 ) , KAB, KABB, KAN, KAT , KATT, KAT 1 , KAT2 , KBI13), 

3 KBA, KBB , K8F , KCH(51), KRL , KTU850), LAST, LB(2), LC ( 3 ) , MQ, 

4 NAJ, NP IT, NRD , NRL , NT, NT J * NURL, NXEP, PTX(99), PTY{99), 

5 RCUI51), RHI3000), RHDOWN, RHUP, RIN, RX, SAU(51>, SEM, 

6 SIZE, U ( 3000 ) , UB i 30001 , URH i 3000) , VA , VAT, VBT , VX(51), 

7 VY { 5 1 ) , XD, XEP , XK, XMP, XQM , XR, XT(1B50),YEP 

C INITIALIZE COMMON AREA TO ZERO 119638 IS THE PRESENT LENGTH ) 
DIMENSION AREMUJ 
DO 7 J=l, 19 638 
7 AREM(JJ=C.C 

C TRANSFER TO DATA INPUT SUBROUTINE 

1 CALL CHN14 

C CHN12 CALCULATES THE POTENTIAL FIELD 

2 CALL CHN12 

C CHN13 CALCULATES TRAJECTORIES AND RHS 

3 CALL CHN13 

C SWITCH TO RESTART FOR NEXT DATA CASE 

6 IF { NURL J 1,2,2 

END 
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SUBROUTINE: MATR I X (N,B , X > 

COMMON AREM, AT X { 5 1 ) , AT Y ( 5 1 ) , AX, AYI51), BASE, CU(51), CUD(5l), 

L DC, DCC , DEL Y , DX, EPS, ER(40), E T X ( 5 I ) , ET Y { 5 1 ) , H , I AS , J AS ( 20 ) , J D 

2 , JOT , JTC3000), KAB, KABB, KAN, KAT , KATT, KATI , KAT 2 , KB(13), 

3 KBA, KBB , KBF , KCHI51), KRL * KTU850), LAST, LB(2), LC(3), MO, 

4 NAJ, NP IT, NRD , NRL , NT, NT J , NURL, NXEP, PTX(99), PTY(99), 

5 RCU ( 51), RH ( 3000 ) , RHDOWN , RHUP, RIN, RX, SAU(5i>, SEM, 

6 SIZE, U ( 3000 ) , UB { 3000 ) , URH i 3000 ) , VA , VAT, VBT , VXC51), 

7 V Y ( 51 ) , XD, XEP , XK, XMP, XQM , XR, XT(1850),YEP 

DIMENSION 6(300) ,X(100) 

M=3*N-2 

B1M+1)=0.0 
B I 2 ) = B ( 2 ) / B 4 1 ) 

X< 1 )=X( 1)/B( 1) 

IF(N-l) 12,12,9 

9 K = 2 

DO 10 , M , 3 

B ( J J =B ( J )-B( J-1)*B( J-2) 

B(J+1)=B(J+1)/BIJ) 

X(K)-(XCK)-B( J-l)*XIK-l> >/B( J) 

10 K = K + 1 

K = K- 1 

DO 11 J= 1 , M , 3 
NB = M— J— 1 
K = K— 1 

11 X(K )=X(K )-B(NB )*X(K+1) 

12 RETURN 
END 
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SUBROUTINE CHN14 
C CHNI4 IS FOR DATA INPUT 

COMMON AREM, ATX l 51 ) , AT Y ( 5 1 ) , AX, AY(51), BASE, CU<51), CULH5L), 

1 DC, DCC, DELY, OX, EPS, ERI40), E TX ( 5 1 ) , ET Y ( 51 ) , H , I AS , JAS ( 20 ) , J D 

2 , JOT , JTI3000), KAB, KABB, KAN, KAT, KATT , KATi, KAT2 , KB(13), 

3 KBA, KBB , K8F , KCHI51), KRL , KTII850), LAST, LB(2), LC(3), MO, 

4 NAJ, NPIT, NRD, NRL , NT, NT J , NURL, NXEP, PTX<99), PTY(99), 

5 RCU (51), RH( 3000 ) , RHDOWN , RHUP , RIN, RX, SAU(5I), SEM, 

6 SIZE, U13000), UB ( 3000) , LRH { 3000) , VA , VAT, VBT , VX15L), 

7 VY( 51 ) , XO, XEP , XK, XMP, XQM , XR , XT(1850),YEP 

READ! 5,100 NP I T , NTP , KAT , KA TT ,NT,NTJ, NTA ,KAN , KBA , KAB , KAT I , KAT2 , IAS 


KEAOI 5, 100 )NPL , KBB 
NT J=NT J- I 

C NPIT =N EGAT I VE , THE PROGRAM STOPS 

C NPIT=0, FIRST GUESS FOR POTENTIAL FIELD IS READ IN FROM GUESS 
C NP I T= + , POTENTIAL FIELD IS READ IN FROM BCREAC 
C NTP=M AX I MUM JT TYPE NUMBER +4 

C KAT = F IRST LINE IN FIRST GRID FOR IMPINGEMENT TEST 
C KAT T=LAST LINE IN FIRST GRID FOR IMPINGEMENT TEST 
C NT=TOTAL NUMBER OF POINTS 
C NTJ = NUM 8ER OF TRAJECTORIES 

C NTA=NUMBER OF INPUT COORDINATES FOR EMITTER 
C KAN=TEST POINT FOR RHS 
C KBA=T EST POINT FOR EQU I POTENT I AL 

C KA8=NUMBER OF LINES TO TRAVERSE TO OBTAIN £ GUIPOTENTI AL LINE FOR CU 
C CALCULATION 

C KAT 1 = F IRST LINE IN SECOND GRID FOR IMPINGEMENT TEST 
C KAT2=LAST LINE IN SECOND GRID FOR IMPINGEMENT TEST 
C I AS=0 FOR 2-0, IA S= 1 FOR 3-0 
C NPL=NUMBER OF POINTS PER LINE 

C K8B=0»— CURRENT DENSITY IS SPACE CHARGE LIMITED 
C KB B=+ CURRENT DENSITIES ARE SPECIFIED TUBEWISE BY RCU REAO 
C NRL=MAX NUMBER OF CYCLES TO CONVERGE ON POISSCN SOLUTION 
NRL = 5 
KRL = NRL 

C KBF=UPPER SYMMETRY SWITCH 
KBF= I 

C RHUP= UPPER BOUND ON KH 
RHUP=0 .0 

C RHDOWN= LOWER BOUND ON RH 
RHD0WN=0.0 

C K 8H=R EAD CONTROL FOR ER 
KBH=KATT-KAT+1 
IFIKAT2J16, 16,17 
17 KBH=KBH+KAT2-KAT 1+ 1 

C PRINT CONTROL FOR LAPLACE TRAJ 
16 J0T=200 

C PRINT CONTROL FOR POISSON TRAJ 
NRD=200 
DO 1 J = 1 , 5 

C READ IN HEADING CARDS 
READ (5,101) 

1 WRITE (6,101) 
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KABB=NTA— 1 
LB I 1 ) =0 
LB ( 2 ) =0 
LC ( 1 )=NT/NPL— 1 
LC ( 2 ) =NPL+ 1 
LC ( 3 ) = NPL 

C JAS IS A VECTOR WHICH DETERMINES TEST POINTS TO BE CHECKED AGAINST 
C EMITTER POTENTIAL 

C JASI 1 )=NUMBER OF POINTS TO BE CHECKED AGAINST EMITTER POTENTIAL 
C JAS(2-( JASt 1)+1))=MESH POINT NUMBERS TO BE CHECKED 
READ ( 5 , 100) ( JA S ( J J » J =1 *20) 

C YEP=PERMITT IVEITY OF FREE SPACE 
YEP = 8 .854E-12 

C RX=RHS SUPRESSION FACTOR 
RX= .4 
EPS= • 1 

C LAST =LOWER SYMMETRY TEST 
LAST = 0 

READ (5,103) AW,VA,H 
C AW = ATOMIC WEIGHT OF ION 
C VA=EM ITTER POTENTIAL 
C H=MESH SIZE 

XQM = S .649E 7 

C XQM=CHARGE TO MASS RATIO OF IONS 
XGM=XQM/AW 

C CC AND DCC ARE USED ONLY TO SUM THE CURRENT AT THE APPROPRIATE STEP 
C DC=DI STANCE FROM 0 TO RIGHT HAND SIDE OF ACCEL ELECTRODE 

C CCC EQUALS CISTANCE BETWEEN ACCEL AND DECEL ELECTRODES 

DC=FLOAT(KATT)*H 
DCC = FLO AT (KAT2-KATT)*H 
IFIKAT2.EQ.0) DCC=C. 

R EAL ( 5 , 1031 VAT, VB T , SIZE 

C VAT = LARGEST POTENTIAL IN EQUI POTENTI AL PRINTOUT 
C VBT = SMALLEST POTENTENT IAL IN E QUI POTE NTI AL PRINTOUT 
C SIZE = STEPSIZE FOR E QUI POTEN T I A L PRINTOUT 
XC=NPL- I 

C RIN = WIDTH OF REGION 
R IN=XC*H 
MO = I 

IFIKBBJ 35,35,36 
36 NN=NTJ+ I 

C RCU= VALUES IF FIXED CURRENT DENSITY 

READ ( 5, 102) (RCU(J) , J = 1 , NN ) 

35 CALL XTCALtNTP) 

C ATX=X- COORDINATES OF EMITTER 

READ (5,103) (ATX(J) , J = l , NTA J 
C AT Y = Y— COORD I NATES OF EMITTER 

READ (5,103)(ATY(J) ,J=l,NTA) 

C ER = T £ S T POINTS FOR IMPINGEMENT 
IF(KAT.EQ.O) GO TO 40 
READ ( 5, 103) I ER ( J) , J=1 ,KBH) 

C CHECK UPPER SYMMETRY 
40 IF(ATY( l)i 2,2,3 

2 KBF=-KBF 

C CHECK LOWER SYMMETRY 

3 IF(ABS(ATY(NTA)-RIN)-1.E-05*H) 30,30,37 
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37 LAST =- 1 

30 DO 4 J=i t NT 

4 J T ( J ) = 0 

C JT IN T SETS UP JT-ARRAY 
34 CALL JTINTINTP) 

C XR = SPECTRAL RADIUS 

12 READ (5,103) XR 

C CGND IT IGNAL EIGENVALUE CALCULATION 

IFU XR*(XR~1. J ) .LT.O. ) GO TO 19 

13 CALL E VC 

C COND ITIUNAL EXIT 

19 IF INPIT) 25,14,27 

C GUESS INITIALIZES THE POTENTIAL FIELD 

14 CALL GUESS 

GO TO 18 

27 CALL BCREADIUI 1 ),UINT)I 

18 XM = LC (31 

XN = LC I 1) 

C XMP^ERROR FACTOR RELATING MESH SIZE 

XMP=. 5*1 XN*XM 1**2/ I XN*XN+XM*XM) 

C NURL=NUMBER OF ITERATIONS ON POISSON SOLUTION 
NURL=3C0 

C KB IS NOW PRINT CONTROL FOR POTENTIAL AND RHS 
KB ( 11 = 1 
DO 15 J = 2 , 13 

15 KB I J ) =0 

DO 22 J = 1 » NT 
UB! J)=U( J ) 

RH ( J ) =0 . 0 

IFI IABSI JT( J)) .LE.8) JT ( J)=G 
IF I J T I J ) ) 3 1 , 22 , 22 
31 IFI LSI 1) ) 21,21,29 

21 LB I 1 ) = J 

29 LB I 2 ) 

22 CONTINUE 

WR I TE ( 6 , 10 4 ) 

RETURN 

25 CALL EXIT 

100 FORMAT ( 1415} 

101 FORMAT ( 7 2H 

1 ) 
102 FORMAT! 7E10. 5) 

103 FORMAT { 7F 10 • 5 ) 

104 FORMAT { 1H0,47X,17H LAPLACE SOLUTION ) 

END 
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C XTCAL CALCULATES COEFF FOR FINITE DIFFERENCE EQN 

SUBROUTINE XTCAL ( N ) 

COMMON AREM, ATX ( 51 ) * AT Y < 5 1 ) * AX* A Y ( 51 ) * BASE, CU(51), CUD(51), 

1 DC, DCC , DELY, DX, EPS, ERI40), E T X ( 5 1 > , ET Y ( 5 1 ) , H , I AS , J AS ( 20 ) , JD 

2 * JOT , JTI30G0), KAB, KABB, KAN, KAT, KATT , KAT1 , KAT 2 , KB(13), 

3 KBA, KBB* KBF * KCHI51), KRL* KTC1850I, LAST, LB(2), LC(3), MO, 

4 NAJ, NP IT, NRD, NRL , NT, NT J, NURL, NXEP, PTX(99), PTY<99), 

5 RCUC51J, RH ( 3000 ) , RHDO WN , RHUP* RIN, RX, S AU ( 5 1 ) , SEM, 

6 SIZE, U ( 3000 } , UB l 30001 , LRH ( 3000 ) , VA , VAT, VBT , VX(51), 

7 V Y { 5 1 ) , XC, XEP, XK, XMP, XQM , XR , XT<1850),YEP 

C XF1S I STANC E UP <-Y DIRECTION) 

C XF I =0 IMPLIES NORMAL DERIVATIVES IN UP DIRECTION 

C X F 1 LESS THAN ZERO IMPLIES THAT THE XFS ARE LOADED INTO THE XT 

C SLOTS DIRECTLY AFTER THE SIGN IS CHANGED ON XF1 
C XF2=CIST ANC E RIGHT (+X DIRECTION) 

C XF 2 = 0 IMPLIES NOR DERIVS IN +X DIRECTION 

C X F 3=0 1ST ANC E OQWN { +Y DIRECTION) 

C XF3 =0 IMPLIES NORMAL DERIVATIVES IN DOWN DIRECTION 

C XF4S 1ST ANC E LEFT <-X DIRECTION) 

C XF4S IMPLIES NOR DERIVS IN -X DIRECTION 
C R = R AD I US FOR AXISYMMETRIC 

C XT ( J ) =COEFF UP 

C XT ( J 4 1 ) SOEFF DOWN 

C XT < J +2 ) =C0 E FF RIGHT 

C XT ( J 4 3 ) =C0 EF F LEFT 

C XT ( J +4 ) SOEFF RH S 

C KTIJ) IS RELATIVE NUMBER FOR XT<J) 

C KT ( J ♦ 1 ) IS THE RELATIVE NUMBER FOR XT(j4l) 

C KT ( J +2 ) IS THE RELATIVE NUMBER FOR XTIJ42) 

C KTIJ43) IS THE RELATIVE NUMBER FOR XTIJ43) 

WR I TE ( 6 , 102) 

IFIIAS) 1,1,2 
C 2- D CALCULATION 

1 DO 10 J = 9 , N » 5 

ISW = 0 
KEL =0 
K = J 4 3 

L = K 4 I 
JNUM = J 

READt 5,100 ) i KT( I ) , I=J ,K) , XF 1 , XF2 , XF3 , XF4 , R 

IF(XFl) 14,15,15 
C STORE COEFF DIRECTLY 
14 XT(J)=-XF1 

XTC J41)=XF2 
XT( J42)=XF3 
XT I J +3 ) =XF 4 
XT ( J 4 4 ) =R 
JNUM=— J NUM 
GO TO 10 
15 XF 1 = XF 1 /H 

XF 2 = XF 2 /H 
XF 3 = XF 3/H 
XF4=XF4/H 

C CHECK FOR NUR DERIVS IN VERTICAL DIRECTION 
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IF ( XF 1 + XF 3— 1 . ) 16,16,17 

16 IF { XF 1 J 18, 18,15 

18 XF 1 = 1 • 

I 1 

GO TO 17 

19 XF3= 1 « 

ISW=-1 

17 XF5 = XF 1 + XF 3 

C CHECK FOR NOR DERIVE IN HORIZONTAL DIRECTION 
IFUXF2+XF4).GT.l.) GO TO 23 
1F< XF2 .EQ.O. ) GO TO 24 
XF4 = 1 • 

KEL=-1 
GO TO 23 
24 XF 2= 1 • 

KEL^l 

23 X F6 = X F 2 + XF 4 

XF7=XF1*XF3+XF2*XF4 

XT { J ) = XF 2* XF 3* XF 4/ (XF5*XF7) 

XT I J + 1 J=XFi*XF2*XF4/4XF5*XF7) 

XT ( J + 2 J = XF l*XF3*XF4/( XF6* XF 7 ) 

XT( J+3)=XFl*XF2*XF3/( XF6*XF7) 

XT{ J+4)=XF1*XF2*XF3*XF4*. 5/XF7 
IF ( I S W ) 20,25,21 
21 XT(J+1)=XJ(J+1)+XT(J) 

XT ( J ) =0. 

GO TO 25 

20 XT(J)=XT(J+1)+XT( J) 

XT ( J+ 1 j =0* 

25 I F i KEL ) 26,10,27 

26 XT I J+2) = XT{ J+2) + XT( J + 3) 

XT( J + 3) = 0.- 
GO TO 10 

2 7 XT< J*3)=XT{ J+3 ) + XT{ J + 2) 

XT( J+2)=0. 

10 WRITE(6, 103 J (KTU ) , I=J ,K) , i XT (I) ,I=J,L) , JNUM 

12 RETURN 

C AXISYM CALCULATION 
2 DO 11 J = 5 , N , 5 

K = J + 3 

JNUM= J 
L = K + 1 
KEL-0 

REAO( 5, ICO) (KTl I ) , I=J ,K) , XF 1 , XF 2 , XF3 , XF4 , R 

IF { XF 1 1 7,9,5 
9 XF 1 = XF 1/H 

XF2=XF2/H 
X F3 = XF 3/ H 
X F 4 = XF A/H 
XF5=XF1+XF3 

C NORMAL DERIVATIVE CHECK 

I F ( (XF2 + XF4J.GT. 1. ) GO TO 28 
IF( XF2 *EQ . 0 • ) GO TO 29 
XF4= 1 « 
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KEL=~i 
GO TO 28 
29 XF 2= 1 • 

KEL = 1 

28 X F6 = XF2 + XF 4 

HR = • 5*R 
hH= • 5*H 
FE = . 1 2 5*H 
IF ( R ) 3 f 3 ♦ 1 3 

C NORMAL DERIVATIVE CHECK 
13 IF{ XF 1 + XF3— 1# ) 4,4,22 

22 XF7=XF1*XF2* XF3*XF4 

XF8 = HR* ( XF2*XF4+XFl*XF3) /XF7+HE*C XF 1 -XF3 )/ (XF2*XF4) 

XF8=XF8*XF5*XF6 

XT ( J )=IR/XF1+HH)*(.5*XF6/XF8) 

XT< J + 1) = <R/XF3-HH)*( . 5*XF6/XF8) 

XT (J+2 )=( HR + HE* ( XF 1— XF 3 ) )*XF5/1 XF2*XF8) 

XT ( J + 3 } = ( HR+HE* ( XF1-XF3) )*XF5/ (XF4*XF8) 

XT { J+4)=«25*R*XF5*XF6/XF8 
GO TO 30 

3 XT ( J + i ) = 0 • 

XL=XF1*XF1*.125 

XF8=( .25+XL/(XF2*XF4) )*XF6 
XT < J 1= XF 6/ XF 8* • 25 
XT ( J +2 ) = XL / l XF 2* XF 81 
XT{ J + 3 )=XL/I XF 4* XF 8 ) 

XT(J+4)=# 5*XF6*XL/XF8 
GO TO 30 

4 IF(XFl) 5,5,6 

5 XF7=R-XF3*H*.25 
XF8=.5*XF6*(R/XF3-HH+XF3/ ( XF 2* XF4 ) * XF 7 ) 

XT ( J ) =C • 

X T ( J * 1 )=.5*XF6*lR-XF3*HH}/< XF3*XF8) 

XT ( J+2)=.5*XF3*XF7/l XF2*XF8) 

XT( J+3)=.5*XF3*XF7/{ XF4*XF8) 

XT! J + 4 )=.25*XF3*XF6*XF7/XF8 
GO TO 30 

6 XT(J+1J=0. 

X F 7 =R +XF i*H* .25 

XF8 = . 5*XF6* i R/XF 1 +HH+ XF 1* XF 7 / ( XF2 *XF4 ) ) 

XT ( J )=.5*XF6*(R+HH*XF1 ) /( XF1*XF8) 

XT( J + 2 ) = .5*XFl*XF7/< XF2*XF8) 

XT ( J+3 )=.5*XF1*XF7/(XF4*XF8) 




XT ( 

J 

+ 4 

) = • 2 

5*XF1* 

XF6* XF 7 /XF 8 

30 

IF! 

[ KEL 

) 

3 

2, 11 

,31 



31 


X T ( 

J 

+ 3 

) = X T 

< J+3>+ 

XT 

(J + 2) 



XT ( 

J 

+ 2 

) = 0 • 





GO 

TO 

1 

1 
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X T ( 

J 

+ 2 

I = X T 

( J+2> + 

XT 

(J + 3) 



XT { 

J 

+ 3 

) = 0 • 





GO 

TO 

I 

1 





7 


XT ( 

J 

) = 

- XF 1 






XT i 

J 

+ 1 

) = XF 

2 





X T ( 

J 

+ 2 

) = XF 

3 





XT i 

J 

+ 3 

) = XF 

4 





XU 

J 

+ 4 

) =R 
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JNUM=— JNUM 

11 WRITE (6, 103) IKK I) .I=J,K> , ( XT( I) ,I=J,L) , JNUM 


GO TO 12 

100 FORMAT 14I5»5F10.0) 

102 FORMAT ( 1 11H KTIJT) KTUT+l) KT(JT+2) KTIJT+3) 

1XTIJT) XT1JT+1) XT1JT+2I XT(JT+3> XTIJT+4) 

103 FORMAT (41 10 » 10X » 5F 10.7*1 10) 

END 
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C JTINT SETS LP THE JT ARRAY 
SUBROUTINE JTINTINTP) 

COMMON AREM, A TX l 51 ) , A TY ( 5 1 ) , AX, AYI51), BASE, CU(51), CUOI51), 

1 DC, DCC , OELY, OX, EPS, ER140), E TX ( 51 ) • ET Y ( 5 1 ) , H • I AS , J AS ( 20 ) , JO 

2 , JOT , JT I 3000 ) , KAB , KABB, KAN, KAT , KATT , KAT1 , KAT2, KB113), 

3 KBA, KBB, KBF , KCH15U, KRL, KT11850), LAST, LB(2), LC<3), MO, 

4 NAJ, NPIT, NRD, NRL , NT, NT J , NURL , NXEP, PTXI99), PTY(99), 

5 RCU (51), RH( 3000) , RHOOWN, RHUP, RIN, RX, SAU(51), SEM, 

6 SIZE, U ( 3000 ) , UBI3000), CRH i 3000 ) , VA , VAT, VBT , VX(51)» 

7 V Y I 5 1 ) , XO, XEP , XK, XMP, XQM , XR, XT(1850),YEP 

KB l 131=0 

IF I IAS.LT.1) GO TO 50 
C AXISYM JT-SETUP 

READ ( 5 , l CC ) JA 
C J A=NUMBER OF CARDS 

C JP=PO INT S PER COLUMN 

JP=LC ( 3) 

JG=NT + 1 
JK= 1 

JD = JK 

DO 11 J= 1 , J A 

READ 15«100)KA»IKB1K) ,K = 1 ,12) 

KC= 1 

7 JE=KB ( KC ) 

IF(JE) 11,11,6 

8 JF = KB ( KC + 1 ) 

C ERROR CHECK ON JT 

JFF= IAB S < JF ) 

1F1JFF.LT.9) GO TO 113 

IFIMODI JFF-4,5 J.NE.O.OR. JFF.GE .NTP) GC TO 111 
113 DO 9 L= 1 , JE 

JT ( JD )=JF 

9 JD = JD + JP 
KC=KC+2 

IF ( JD- JG J 7,12,12 
12 JK = J K + 1 

JD = JK 
GO TO 7 

11 CONTINUE 

READ( 5, 10DBASE 
31 kETURN 

C 2-0 JT SETUP 
50 READ 15,100) JA 
J B = 0 

60 DO 1 10 J= 1 , J A 

READ I 5,10C)KA»IKB(K) ,K=1 ,12) 

KC = 1 

DO 10 K— 1, KA 
70 JD = KB ( KC ) 

IF ( JD ) 10,10,60 

80 JE = KBiKC+1) 

C ERROR CHECK ON JT 
JFF= IAB S ( JE ) 
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IFUFF.LT. 9) GO TO 114 

IFiMOOi JFF-4, 5) .NE . C. OR . JFF . GE . NTP ) GO TO 111 
DO SOL = 1 # JO 
JB = JB+1 
C JT=POINT-TVPE VECTOR 
SO JTIJB) = JE 

KC * KC 42 
GO TO 70 
10 KC * 1 

110 CONTINUE 
GO TO 31 

111 WR I TE < 6 » 1021 
WRITE(6tl00)KAt(KBCK) »K = 1 ,131 

CALL EXIT 

100 FORMA T I 1415) 

101 FORMAT ( F 10. 5 ) 

102 FORMAT ( 18H ERROR IN JT DATA. ) 

END 


38 



u u u o 


C EVC CALCULATES THE SPECTRAL RADIUS UF THE ITERATION MATRIX 

SUBROUTINE EVC 

CCMMON AREM, ATXl 51 ) ,ATY< 51) , AX, AYI51), BASE, CU(51>, CUDI51), 

1 DC, DCC , DELY, DX, EPS, ER140), ETX I 5 1 ) , ET Y I 51 ) , H , I AS , JAS I 20 > , JD 

2 , JOT , JTI3C00), KAB, KABB, KAN, KAT, KATT , K AT 1 , KAT2 , KBI13), 

3 KBA, KBB , KBF, KCHI51), KRL , KTI1850), LAST, LB ( 2 ) , LC ( 3 ) , MO, 

4 NAJ, NP IT, NRD, NRL , NT, NT J , NURL , NXEP, PTX199), PTY(99), 

5 RCUI51), RHI3000), RHDOWN , RHUP, RIN, RX, SAUI51), SEM, 

6 SIZE, UI3000), UBI3000), CRH ( 3000 ), VA , VAT, VBT , VXI51), 

7 VY ( 5 1 > , XD, XEP , XK, XMP, XQM , XR, XTI1850),YEP 
DIMENSION AA I 300 ) , X( 100) 

EQUIVALENCE (URHI 1 ) , XI l) ) ,(URH(101) ,AA(1)) 

THE ABSOLUTE VALUE OF THE LARGEST EIGENVALUE OF THE MATRIX IS ITS 
SPECTRAL RADIUS . THE TERMS EIGENVALUE AND SPECTRAL RADIUS 
WILL BE USED INTERCHANGEABLY 
INITIALIZE VECTOR ITERATES 
DO 3 J = 1 , NT 
RHI J ) =0. C 
3 U(J)=C.O 

JS=- 1 

C INITIALIZE EIGENVECTOR 

DO 8 J D= 1 , N T 
IFIJT(JD)-8) 6,6,7 

6 U I JD ) -C » 

GO TO 8 

7 U ( JD ) = 1 . 

8 CONTINUE 

C NUMBER OF COLUMNS 
NL I N=LC I 1) + 1 
J M = LC ( 3 )— 1 

C BEGIN MINMAX PROCESS 

9 DO 22 KK=1 , ICOO 

C CALCULATE NEW VECTOR 

DO 29 ML = 1 , NL I N 
C FIRST MESF POINT 

KZ=IML-1)*LCI 3 ) +1 
C LAST MESH POINT 
KC = KZ + JM 
J V= 1 
JU=1 

DO 38 K = KZ , KC 
JZ = J T ( K ) 

KU=KT< JZ)+K 
KD=KT( JZ+l )+K 
KR=KTI JZ*2UK 
KL = KTI JZ + 3) + K 
IFIJZ.GT.8) GO TO 37 
AAI JU)=0. 

AA( JU+1)=1. 

AAI JU+2)=0. 

SUM=0. 

GO TO 35 

3 7 SUM = U(KR)*XT (JZ+2)+-UIKL)*XT IJZ + 3) 

AAI JU)*-XT( JZ) 

AAI JU4l)=l. 
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AA I J U + 2 ) =-XT (JZ+1) 

I FI (K.EC.KZ).OR. U JTIK-l) .GT.O) .AND. ( KT < J Z ) . EQ . -1 ) ) ) GO TO 36 
AAl JU)=0. 

36 IF I (K.EC.KC J .OR. ( { JT1K+1) .GT.O) .AND. I KT I JZ+1) .EC.l) ) )GO TO 35 

AA ( JU+2 ) =0. 

35 X( J V ) = SUM 

JV=JV+1 
38 J U= JU + 3 

CALL MA TR I X ( LC I 3 3 * AA ( 2 ) * X ) 

J V= 1 

CO 3 A K =KZ » KC 
RH( K J = X1 JV 3 

IFUTIK3.LE.fi) RH IK) =0. 

34 JV= J V+ 1 

29 CUNT INGE 

C ThE MATRIX MOST BE APPLIED TWICE AS ITS TWO-CYCLIC NATURE 
C REORDERS THE MESH POINTS AFTER ONE MULTIPLICATION. 

C THE SPECTRAL RADIUS SQUARED IS THE RESULT OF THE MINMAX 
C PROCEDURE AND THUS ITS SQUARE ROOT MUST BE TAKEN 

C SWITCH ALLOWING MATRIX TO BE APPLIED TWICE 

IFUS) 13,15,15 
13 DO 14 K= 1 ,NT 

C STORE OLD VECTOR IN UB (J ) 

UB ( K ) =U( K ) 

14 U(K)=RHJK) 

GO TO 22 

C INITIALIZE FOR RATIO TEST 

15 XL-0 . 0 
XS=1.0 

DO 20 JD= 1, NT 
IFIUI JD))2C,20,16 
16 X=RH(JD)/UB(JD) 

IF(XL-X) 17,18,18 
C XL IS MAXIMUM RATIO 

17 XL=X 
NNL=JD 

18 IF(XS-X) 20,20,19 

C XS IS MINIMUM RATIO 

19 XS=X 
NS = JD 

20 CONTINUE 
YL=RHI NNL ) 

DO 21 JD= 1 * NT 

C SCALE SO THAT VECTOR DOESNT GET TOO LARGE IN MAGNITUDE 

21 U1 JD)=RH( JD 3/YL 

C CONVERGENCE TEST ( IE-4 IS ARBITRARY) 

I F ( X L- X S— 1 • OE— 4 ) 24,24,22 

22 J S = — J S 

23 WRITE (6,102)XS,XL 

24 XR=SQRTI.5*(XS+XL1) 

25 WRITE(6, 1C1 ) XR , KK 

C ERROR CHECK 

IF ( (XR.GE .1. ) .OR. (XR.LE.O. ) ) GO TO 28 

27 RETURN 

28 CALL EXIT 

100 FORMAT! 1415) 
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1C1 FORMAT {4HOXR=F11.8,2H I5.39H ITERATIONS REQUIRED TO CONVERGE 0 

*N XR ) 

102 FORMAT ( 7H0XS XL 2F15.8) 

END 
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C SUBROUTINE GUESS GENERATES A FIRST GUESS 
SUBROUTINE GUESS 

COMMON AREM, ATXI 51 ) ,ATY ( 5 1 ) , AX, AY151), BASE, CUI51), CUDI51), 

1 DC, DCC , DELY, DX, EPS. ERI40), ET X ( 51 ) , ETY ( 51 ) , H , I AS , J AS ( 20 ) , J D 

2 *JOT» JTI3C00), KAB, KABB, KAN, KAT, KATT , KAT1 , KAT2, KB(13)» 

3 KBA, KBB, KBF , KCHI51), KRL , KK1850I, LAST, LB(2), LC ( 3 ) , MO, 

A NAJ, NP IT, NRD, NRL , NT, NT J , NURL, NXEP, PTX199), PTYI99), 

5 RCUI51), RH ( 3000 ) , RHDOWN , RHUP, RIN, RX, SAU(51), SEM, 

6 SIZE, U13C00), UB (3000) , LRH t 3000 ) , VA , VAT, VBT , VX(5i), 

7 V Y < 5 1 ) » XD, XEP , XK, XMP, XQM , XR, XT(1850),YEP 
DIMENSION GEP t 7 ) 

EQUIVALENCE i UB t 1 ) ,GEP ( 1 ) ) 

C GEP ( 1 ) IS ASSOCIATED WITH JT-TYPE 0,GEP(2) WITH TYPE1 , ETC. 

READ! 5, IOC) (GEPI J ) , J=1 , 7) 

DO 1 J= 1 »NT 
U( J ) =0 . c 

IF( IABSIJTIJ ))-9) 2,1,1 
2 K=JTIJ)+1 

U( J )=GEP ( K ) 

JT ( J ) =0 

1 CONTINUE 

RETURN 

100 FORMAT ( 7F10. 5) 

END 
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SUBROUTINE CHNI2 

C CHNI2 CALCULATES THE POTENTIAL FIELD 

COMMON AREM, ATX( 51 ) ,ATY( 51 ) , AX, AY(51), BASE, CU151), CUDI51), 

1 DC, DCC, DELY, DX, EPS, ER(4C), ETX l 5 1) , ET Y ( 51 ) , H , I AS , J AS ( 20 ) , JD 

2 , JO T , JT ( 3000 ) , KAB, KABB, KAN, KAT, KATT , KAT 1 , KAT 2 , KB(13), 

3 KBA, KBB, KBF , KCH151), KRL , KK1850), LAST, LB<2), LC ( 3 > , MO, 

4 NAJ, NP IT, NRD, NRL, NT, NT J , NURL, NXEP, PTXI99), PTY(99), 

5 RCU ( 5 1 ) , RHI300C), RHDOWN , RHUP , RIN, RX, SAU(51), SEM, 

6 SIZE, U ( 3C00 ) , UB I 3000 ) , URH (3000) , VA , VAT, VBT , VX(51), 

7 V Y ( 51) , XD, XEP , XK, XMP, XQM , XR , XT(1850),YEP 

DIMENSION AA( 300) ,X(100) 

EQUIVALENCE (URH ( 1 ) ,X( 1) ) , ( URH (101 ) ,AA( 1) ) 

SX=1.0 

XM=.25*XR**2 
NL I N=LC ( 1) + 1 
JA = LC ( 3 i — 1 
1 XW= 1 . 

DO 23 NN = 1 * NURL 
XEP =0 . 

ML= 1 

2 DO 29 LL =ML , NL I N , 2 
KZ = ( LL- 1 ) *LC ( 3 ) + 1 

KC=KZ+ JA 
J V= 1 
JU= 1 

00 38 K = KZ , KC 
JZ = J T ( K ) 

KU= K + K T ( JZ ) 

KD = K +K T ( JZ + 1 ) 

KR=K+KT ( JZ+2 ) 

KL=K+KT ( JZ+ 3 ) 

IF(JZ.GT.O) GO TO 37 
AA ( J U ) =0 . 

AA ( JU ♦ 1 ) = 1 . 

AA ( J U + 2 ) =0 . 

SUM = U ( K ) 

GO TO 35 

3 7 SUM = U ( KR )*XT(JZ+2)+U<KL)*XT(JZ+3>+RH(K)*XT ( JZ+4 ) 

AA(JO)=-XT(JZ) 

AA( JUtl )=1. 

AA(JU+2)=-XT(JZ+I) 

IF( CK.EQ.KZ) .OR. I( JT(K-1 ) .GT.O) .AND. (KT( JZ) .EQ.-l ) ) ) GO TO 36 
AA( JU)=0. 

SUM=SUM+U(KU)*XT ( JZ) 

36 IF( ( K .EQ .KC ) .OR. ( ( JT< K+l ) .GT.O) .AND. ( KU JZ+1) .EQ.l) ) ) GO TO 35 

AA ( J U + 2 ) =0 . 

SUM= SUM +L(KD)*XT( JZ+1) 

35 X( J V ) =SUM 

JV=JV+1 
38 JU= JU+ 3 

CALL MATRIX(LC(3),AA(2) ,X) 

J V= 1 

DO 34 K =KZ , KC 
DIF=X(JV)-U(K) 

U(K ) = XW*D IF +U ( K ) 
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DIF =ABS(DIF } 
IF(DIF.LE.XEP) GO TO 34 
XEP=DIF 
NXEP=K 

34 JV=JV+1 

29 CONTINUE 

IFINN.GJ.l) GO TO II 
XW=1./I 1.-2. *XM) 

GO TO 12 

11 XW=1 ./I l.-XM*XW) 

12 IFIML.GT.l) GO TO 3 

ML= 2 


GO TO 2 

3 IF(LBll)) 18,18,15 

15 KG=LB ( 1 ) 

C CALCULATION OF NEGATIVE JT-TYPE MESH POINTS 
KH=LB! 2 > 

00 17 JD=KG» KH 
KE=— JT I JO ) 

IF I KE ) 17,17,16 

C CONTRIBUTING POINTS ARE DETERMINED 
16 KU = K T ( KE I + JD 

KO=KT ( KE+1 J *JD 


KL=KT ( KE+2) +JD 
KR = KT ( KE + 3 } +JD 

UI JD ) =X T ( KE *4 ) *RH ( JD ) +X T ( KE ) *U I KU ) + XT ( KE «- 1 I *U ( KD ) +X T ( KE + 2 ) *U ( KL 
1 )+XT(KE*3)*U!KRJ 
17 CONTINUE 

C CONVERGENCE TEST 
18 IFIXEP*XMP.LE.EPS) GO TO 24 

23 CONTINUE 

24 KNUM=JAS(1) 

IFIKNUM.LE.O) GO TO 27 

C CHECK IF ANY POTENTIALS ARE ABOVE THOSE OF THE EMITTER 
DO 28 J=2,KNUM 
JN= JAS ( J ) 

C MULTIPLICATION BY XQM IS FOR SIGN CONVERSION ONLY 
IF! (UI JN )-VA)*XQM) 28,28,25 
28 CONTINUE 

GO TO 27 

25 SX=SX#RX 

C KAN IS THE TEST POINT 

WRITE 16, 101 ) JN,U( JN) ,SX 
DO 26 J = 1 * NT 
UI J ) = UB I J 3 

26 RH( J )=RHI J )*RX 
RHUP=RHI KAN 3 

GO TO 1 

27 CALL TEST(NN) 

RETURN 

101 FORMAT! 3 8HOPO I N T / VALUE /TOTA L SUPPRESSION OF RHS I5.2E15.6) 

END 
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C PRINTOUT OF POTENTIAL FIELD 
SUBROUTINE TWOUT(KK) 

COMMON AREM, ATX( 51 ) , ATY( 51 ) , AX, AY(51), BASE, CU(51), CUD151), 

1 DC, CCC, DELY, DX, EPS, ER( 40 , E TX 1 5 I ) , ET Y ( 5 1 ) , H , I AS , J AS ( 20 ) , J D 

2 , JOT , JTI3000), KAB, KABB, KAN, KAT, KATT , KAT 1 , KAT2 , KBII3), 

3 KBA, KBB, K8F , KCH (51), KRL, KTU850), LAST, LB<2), LC ( 3 ) , MO, 

4 NAJ, NPIT, NRD, NRL , NT, NT J , NURL, NXEP, PTXI99), PTY(99), 

5 RCUI51), RH ( 3000 ) , RHDOWN » RHUP, RIN, RX, SAU(51), SEM, 

6 SIZE, UI3000J, UB ( 3000 > , URH C 3000) , VA , VAT, VBT , VX(5l), 

7 VY { 5 1 ) » XD, XEP , XK, XMP, XQM , XR , XT1I850),YEP 

NN=KK 

IF(MO) 2,10,2 

2 WRITE! 6, 101) NN, XEP, NXEP 

13 WRITE (6,102) NT 

WRITE<6, 100) 

J = 1 

DO 4 K=1»NT 
KT(J) =K 

XT I J l=U ( K ) 

J = J«-1 

1 F ( J-9 3 4,3,3 

3 WRITE (6,103)(KT(M) ,M = 1 ,8) , (XT (M) ,M = 1 ,8) 

J = 1 

4 CONTINUE 

5 IFIJ-2) 8,6,6 

6 DO 7 K = J » 8 
KT 1 K ) =0 

7 XT ( K ) =C . C 

WRITE (6,103)(KT(M),M=1,8) , ( XT ( M) »M = 1 ,8) 

8 RETURN 

10 WRITE (6,104) 

C POTENTIAL FIELD IS AVERAGED FOR FINAL VALUES 
JOT =NRD 

c set switch to end problem 

NURL=- 7 7 
DO 11 J= 1 ,NT 

11 U(J ) = .5*(U(J ) +UB { J ) ) 

GO TO 2 

100 FORMAT ( 1H0 , 10X ,20H MESH PCI NT NUMBERS 45X.L7H POTENTIAL VALUES) 

101 FORMAT! 7H0AFTER I5.44H ITERATIONS ON U THE MAXIMUM CHANGE IN U 
♦IS F10.5.32H VOLTS AND OCCURS AT MESH POINT 15) 

102 FORMAT ( 1HCI5,9H U VALUES) 

103 FORMAT ( 1H 8I5,8F11.4) 

104 FORMAT ( 4 7HC U-FIELD IS AVERAGE OF THIS AND PREVIOUS CYCLE ) 

END 
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C TEST CHECKS CYCLES AND PRINTOUTS 
SUBROUTINE TESHKK) 

. COMMON AREM, ATX( 51 ) « ATYI 51 ) « AX, AY(5I), BASE, CU(51), CUDI51), 

1 CC, DCC , DELY, DX f EPS, ERI4C), E TX { 5 1 ), ET Y { 5 1 ), H , I AS , J AS < 20 ), J D 

2 , JOT , JTI3C00), KAB, KABB , KAN, KAT, KATT , KAT 1 , KAT 2 , KB(13), 

3 KBA, KBB , KBF , KCHI51), KRL , KTU850), LAST, LB(2), LC ( 3 ) , MO, 

4 NAJ, NPIT, NRD , NRL , NT, NT J , NURL, NXEP, PTX(99), PTY{99), 

5 RCU (51) , RH { 3000) , RHDOWN, RHUP, RIN, RX , SAU(51), SEM, 

6 SIZE, U ( 3 COO ) , UB (3000) , IRH { 3 000 ) , V A , VAT, VBT , VX(51), 

7 V Y { 5 1 ) , XD, XEP , XK, XMP, XQM , XR, XT(1850),YEP 

NN = KK 

I E { KRL— NRL ) 2,1,2 
1 CALL BCDUMP { U(l) ,U(NT) ) 

NURL = 6 C 

2 IWRL=NRL-KRL+I ABS(MO) 

IF ( KB ( IWRL )) 4,4,3 

C CONDITIONAL PRINTOUT OF POTENTIAL FIELD AND EQU I POTE NT I A LS 

3 CALL TWOIT(NN) 

CALL ECLINE 

C SWITCH FOR LAST CYCLE 

4 I F ( NURL ) 6,5,5 
C TEST ON CYCLES 

5 IF ( KRL ) 7,6,6 

6 RETURN 

7 WRITE (6,101) 

C TRANSFER FOR NEW SET OF DATA 

8 NURL =- 7 7 
RETURN 

101 FORMAT ( 11H1NEXT CASE. 

END 


46 



U O' 


C EGLINE CALCULATES THE EGU I POTE N TI AL S 
SUBROUTINE EGLINE 

COMMON AREM, ATX < 5 l ) , AT Y I 51 ) . AX, AY151), BASE, CU(51), CUD ( 5 1 ) , 

1 DC, DCC , DEL Y , DX, EPS, ERI40), E TX I 51 ) , ET Y ( 5 1 ) , H , I AS , J AS ( 20 ) , J D 

2 , JOT , J T ( 3 COO ) , KAB, KABB , KAN, KAT, KATT, K AT I , KAT2, KB(I3), 

3 KBA, KBB, KBF , KCH(51), KRL , KTU850), LAST, LB(2), LC < 3 ) , MO, 

4 NAJ, NP I T , NRD , NRL , NT, NT J , NURL, NXEP, PTX(99), PTYC99), 

5 RCU (51), RH { 3000 ) , RHDOWN, RHUP, RIN, RX, SAUI51), SEM, 

6 SIZE, U ( 3C00 1 , UB { 3000 ) , LRH i 3000 ) , V A , VAT, VBT , VX<51), 

7 VY(51), XD, XEP , XK, XMP, XQM , XR, XT(1850),YEP 
SUM=FLOAT(LC( I)j*H 

XK = I • 

C NO EQUIPGTENTIALS IF SIZE = 0 
IF I SIZE ) 1,27,1 
1 POT BN = VA T 

WRIT£(6, 1021 
JC=LC ( 1 ) 

JD=LC<3 J-l 
DX = H 

31 I=DX*100.0 

C SCALE MESH SIZE IF NECESSARY FOR PRINTING PURPOSES 
IF 1 I ) 29, 29, 3C 

29 0X=DX*1C.0 
S UM = SUM* 10 • 

R IN=R I N* 10 • 

XK= XK* 10 • 

GO TO 31 

30 IFH R IN.LE . 100. > -AND. ( SUM.LE . 100. ) ) GO TO 41 
XK=XK / 10 . 

DX = DX/1C . 

SUM= SUM / 10 . 

R IN = R IN/ 10. 

GO TO 30 
A 1 K = XK + • 5 

IF(K.EG.l) GO TO 2 
WRITE (6,101 JH,DX,XK 

2 J E = LC I 2 ) 

J ED = J E * J D 
L= 1 

3 AX = 0 • 0 

33 CO 22 J J= L , JC 

K S= 1 

5 AA Y = 0 . 0 

CO 19 K = JE , J ED 

6 IF(KS) 8 * B , 7 

7 M= 1 

J = K— LC ( 3 ) 

GO TO 9 

8 J = K — 1 

CHECK TO SEE IF POTENTIAL IS BETWEEN U(J) AND U(K) 

IF ( (UI K 1-PQTEN) *< U( J J-POTEN) ) 10,10,18 

10 C I F=AB S { 0 ( J )-U(K) ) 

IF ( D IF ) 13,13,11 

11 IF ( M ) 12*14.12 

C LINEAR I NP ERPOLAT ION IN HORIZONTAL DIRECTION 
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VX(L }=A8S(UI Jl-POTEN) /DIF*DX*AX 
VY{ LI^AAY 
GO TO 15 

13 VX 1 L ) =AX +DX 
V Y ( L ) =AA Y 
GO TO 15 

14 VX { L ) =A X +DX 

C LINEAR INTERPOLATION IN VERTICAL DIRECTION 
V Y { L ) = A6S(U( J )-PCTEN)/DIP*DX+AAY 

15 IFCL-6) 17,16,16 

16 WRITE (6,100)P0TENt < VXU ) tVYlI ) ,1 = 1,6) 

L^O 

17 L=L + 1 

18 AA Y=AA Y+DX 

19 CONTINUE 
IF(KS) 21,21, 2C 

20 KS =0 
M=0 

JE=JE+1 
GO TO 5 

21 JE=JE4JD 
JED=JE+JD 
AX=AX+DX 

22 CONTINUE 
IFIL-2) 25,23,23 

23 DO 24 J = L t 6 
VX( J)=C.O 

24 VY ( J ) =0 • 0 

WRITE (6,100) POTEN » ( VX < l ) , V Y ( 2 ) ,1 = 1,6) 

C STEP POTENTIAL DOWN 

25 POTEN=POTEN-SI IE 

IF! PGTfcN-VBT) 27,2,2 
27 R1N=RIN/XK 
RETURN 

100 FORMAT { 17H0 POTENT I AL (X,Y) F8.1,2H 6(2H ( F6 . 3 , 1H , F6 * 3 , 2H ) )) 

101 FORMAT ( 3H H=F11*6,5H DX=F11.6,15H SCALE FACTOR=F 1 1 .6 ) 

102 FORMAT ( 1H045X23H EQUI POTENT I AL PRINTOUT ) 

END 
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SUBROUTINE CHN13 

C CHNI3 CALCULATES THE TRAJ AND RHS FOR 2-0 CASES 

COMMON AREM, ATX( 51 ) ,ATY 1 51 ) , AX, AY{51), BASE, CU(5I), CUD(51), 

1 DC, DCC , DELY, DX, EPS, ER(4C), E TX ( 5 I ) , ET Y ( 51 ) , H , I AS , J AS ( 20 ) , J D 

2 , JOT , JT { 3000 ) , KAB , KABB, KAN, KAT , KATT , KAT1 , KAT 2 , KB(13), 

3 KBA, KBB , KBF , KCH<513, KRL , KT<1850), LAST, LB(2), LC ( 3 ) , MO, 

4 NAJ, NPIT, NRD, NRL , NT, NT J , NURL, NXEP, PTX<99), PTY(99), 

5 RCU151), RH { 3000 ) , RHDUWN, RHUP , RIN, RX, SAU(51), SEM, 

6 SIZE, U I 3000 ) , UB { 3000 ) , URH (3000 ), VA , VAT, VBT , VX<51), 

7 VY (51), XD, XEP , XK, XMP, XQM , XR , XT(1850),YEP 

IF(NRL.EQ.KRL) NPIT=C 

NPIT =NP I T + l 

wR I T E ( 6 , 100 ) NP I T 

SUMTWO=0 .0 

SUMTR I =C . 0 

C STORE LAST ESTIMATES 
DO 35 J=1 , NT 
UB ( J 3 = U 1 J ) 

35 URH ( J ) =RH ( J ) 

C INITIALIZATION 

DO 1 J= 1 , 5 1 

C SAU WILL HAVE EMITTER CURRENTS 
SAU ( J ) =0 . 0 

C ETX.ETY WILL HAVE COORD OF EMITTER SURFACE AT EQUAL ARC INCREMENTS 
ETX ( J ) = 0.0 
ETY(J3=C.O 

C KCh=TRAJECTCRY REFLECTION COUNTER 
KCH1 J 3 =- 1 

C AYCJ) IS Y-COCRD OF J-TH TRAJ 
A Y ( J ) =0 • 0 

C VY(J) IS Y— CCMPONEN T OF VELOCITY OF J-TH TRAJ 
VYI J 3=0.0 

C VX(J) IS X-COMPONENT OF VELOCITY OF J-TH TRAJ 
VX ( J ) = C. 00001 

C CU1J) IS CURRENT IN J-TH TUBE 
1 CU(J)=0.0 

DO 11C0 J= 1 * S9 

PTX.PTY ARE COORD OF E CU I PO TE N T I AL LINE USED TO CALCULATE 
CURRENT DENSITY AT THE EMITTER 
PT Y ( J ) = C.O 
1100 PTX ( J ) =0 ,0 

NJT =NTJ 

C PEQ CALCULATES TEE EQUI POTENT I AL LINE IN ORDER THAT THE CURRENT 
C DENSITIES CAN Bb CALCULATED 
55 CALL PEC 

C ARC DIVIDES THE ARC LENGTH OF THE EMITTER INTO EQUAL INCREMENTS 
C FOR THE CURRENT DENSITY CALCULATION 
54 CALL ARC 

C KBA IS THE POINT NUMBER WHICH WILL BE TAKEN FOR THE EQU I POT ENT I AL 
C LINE 

PO TEN = U ( KBA ) 

DELU=ABS ( VA-PQTEN) 

C X K = CGN ST ANT IN CHILDS LAW FORMULA* Db LTA U**3/2 
XQ =AB S ( XQM ) 

XK=4.C/9.0*YEP*SQRT(2.0*XQ *OELU)*D£LU 
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NAJ=NTJ+1 

IFI JOT . L E • 0 ) GO TO 2 
WRITEIfc, 105) 

C INITIALIZE 
2 AX = 0 »0 

NOT = JO T 
DX = H 
DELY=DX 

C JC=NUMBER OF LINES TO TRAVERSE 
JC= LC I 1 ) 

C J E = F IR ST ME S F POINT IN SECOND LINE 
J E = LC { 2) 

C JC=NUM!B ER OF MESH POINTS PER COLUMN 
J D=LC { 3 ) 

S IGN=XCM/XG 
RM=DX/YEP*SIGN 

IF { IAS.GT.O) RM=DX*RM 
XD= .5*XCM/DX 

6 DO 22 JN= 1 , JC 
JED=JE+JD-1 

7 AX = AX +DX 

C TRCU CALCULATES THE CURRENT DENSITIES AND CALLS TRAJ WHICH 
C CALCULATES THE TRAJECTORIES 
IFI IAS. EG. 0) GO TO 34 
CALL ATRCUI JE , JED) 

GO TO 10 

34 CALL TRCU(JEfJED) 

C CORRCT CONDITIONALLY SHIFTS OR TERMINATES THE TRAJECTORIES 
C AT THE GRIDS 
10 CALL CORRCT I JN ) 

11 IFI IAX-3.0*DX-DC)*IAX~4.0*DX-DC> ) 12,12,14 
C SUMTWO IS THE CURRENT AFTER THE FIRST GRID 

12 N A J=NA J 
SUMTWG=C.5*(CUC 1)+CU(NAJ) ) 

IFI IAS.GT.O) SUMTWO =2. *SUMTWO 
DO 12 K = 3, NA J 

13 SUMTWO^SLMTWO+CUlK-1 ) 

GO TO 17 

14 IFI ( AX-3.0*DX-DC-DCC)*(AX-4.0*DX-DC-DCC) ) 15,15,17 

C SUMTRI IS THE CURRENT AFTER THE SECOND GRID 

15 NAJ^NAJ 

SUMTRI=C.5*I CUI 1 J+CUINA J) i 
IFI IAS.GT.O) SUMTRI=2. ^SUMTRI 
DO 16 K=3,NAJ 

16 SUMTR l = SUMTRI+CU(K-l) 

17 IF(NOT) 19, IS, 18 

18 NOT =N0 T— 1 
NT J=N T J 

WRITE! 6, 101) AX, (K,KCH(K) , A Y I K) , V X { K ) , V Y ( K ) , K= 1 , NT J ) 

IS CO 20 J= JE , JED 

20 RH I J ) =C • 0 

IFI IAS. EG. 0) GO TO 48 
C ACALR CALCULATES RHS FOR AXISYM 
CALL ACALRI JE, JED) 

GO TO AS 
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C C ALR CALCULATES ThE RHS FOR 2-0 

48 CALL CAL R ( J E , JED ) 

49 DO 21 K = JE , J ED 

21 RH ( K ) =R H ( K ) *RM 

22 JE=JE4JD 
IW=NRL-KRL+1 

I F I KBB ) 207,207,208 
2G8 KK = NT J-N J T 

NN= N J T + 1 
DO 206 J = 1 , NN 
J J= J +KK 

206 RCU t J ) = RCU ( J J ) 

207 NT J =N J T 

JOT = I A 6 S ( KB i IW + 1) ) 

IF ( KB ( IW)) 30, 30,29 
29 S UM = 0 . C 

SU = 0 • c 

C THRUST AND POWER CALCULATION 
C U- . 5 $ CU 
NA J=NA J 

IF(IAS.EO.O) CU(NAJ) = . 5*CU( NAJ) 
DO 198 J = 1 , NAJ 
CU( J)=CU( Jd*SIGN 
CUCt J)=CUD( J)*SIGN 
SAUtJ) = S A U I J )* S IGN 

I F ( J-l) 200, 200, 201 

200 V=VX(J) 

GO TO 199 

201 IF(J-NAJ) 20 2, 2C3 , 2 03 

203 V=VX(J-1) 

GO TO 199 

202 V=.5*IVX(J)+VXIJ-1) ) 

199 AY { J I = CUIJ)*V / XQM 

V Y ( J)=.5*V*AY( J) 

SUM= SUM *AY i J ) 

198 SU = SU + V V ( J ) 

IF! IAS.GT.O ) GO TO 41 
WRITE (6,109)(J,AY(J) ,J = 1,NAJ) 
WRITE <6, 110 )SUM,SU 
GO TO 4 C 

41 WRITE(6,116) ( J , A Y ( J ) , J = 1 , NA J ) 

WRI TE< 6, 106) SUM, SU 
40 I F { IAS.GT.O) GO TO 47 

SAU(NAJ)=.5*SAU(NAJ) 

SAUt 1 )=.5*SAU( 1) 

C SUMONE IS THE INITIAL CURRENT 
47 SUMUN E = 0 #0 

DO 8 J = 1 » N A J 

8 SUMONE=SUMUNE+SAU( J) 

WRITE (6,111 J{J,CUO(J) ,J=1, NAJ) 

I F ( IAS.GT.O ) GO TO 39 

WRITE (6, 108 )( J, SAUtJ) ,J = 1 ,NAJ) 

WRITE ( 6, 102 )SUMONE 

ER I = SUMONE/ SEM 

WRITE (6,107 )ERI 

ER I = ER I* SEM 
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GO TO 38 

wR I TE ( £ , 115) (JtSAU(J) » J=1 , NA J ) 

WRITE16, 112) SUMONE 
ER I=SUMONE 
38 ERR =SUM TWO/ SUMONE* SIGN 

IFIDC.LE.O.) GO TO 30 
ERA=ERR*ERI 
ERR=ERR*100.0 
IFC IAS.GT.O) GO TO 45 
WRITE ( 6 * 103 )ERA , ERR 
I F ( CCC . EQ .0 . ) GO TO 30 
ERR = SUMTR I / SUMONE *S I GN 
ERA=ERR*ER I 
ERR = ERR* 100.0 
WRITE ( 6 » 1 04 )ERA , ERR 
GO TO 30 

45 WR I TE ( 6 * 113) ERA, ERR 

IF! DCC .EQ.O. ) GO TO 30 
ERR=SUMTRI/SUMONE*SIGN 
ERA 3 ERR*ER I 
ERR =ERR* 100 . 

WR I T E ! 6 » 114i ERA, ERR 

30 IF(NURL) 32,31,31 

C RTEST CHECKS THE BOUNDS ON RHS 

31 CALL RTEST 

32 RETURN 

100 FORMAT ! 8HCCYCLE 12) 

101 FORMAT I 1HCF8.4, 4X1 2 , 1 3 XI 2 , 9XF 9.4 , 2E1 3. 5/ ( 1 3 X 1 2 , 13X 1 2 , 9X F9 . 4 , 

* 2E13.5I) 

1C2 FORMAT! 23H0T0TAL INITIAL CURRENT 3 E12.6.19H AMPS/(UNIT H) ) 

103 FORMAT < 36H0TRANSMI TTED CURRENT AT ACCEL. GRI D= E12.6,24H AMPS/ 

* ( UN I T H) WHICH IS F6.2.32H PERCENT OF THE INITIAL CURRENT. ) 

104 FORMAT (36H0TRANSMITTED CURRENT AT DECEL. GRID 3 E12.6.24H AMPS/ 

* I UN IT HI WHICH IS F6.2.32H PERCENT OF THE INITIAL CURRENT. ) 

105 FORMAT (47H0X-C00RD TRAJ NUM REFLECTION COUNTER Y-COORD 

* 3X, 23HX— VEL COMP Y-VEL COMP ) 

106 F0RMATI24H0T0TAL THRUST (NEWTCNS) E14.6//21H TOTAL POWER (WATT 

*S) E14.6) 

107 FORMAT! 33H0AVER AGE INITIAL CURRENT DENSITY 3 £12.6, 17H AMPS/fUNI 

*T H)**2 I 

108 F0RMATI32H0INI TIAL CURRENTS (AMPS/UNIT H) / / C 7 ( 1H I2,E14.6))) 

1C9 FORMAT (53HOTHRUST DISTRIBUTION BY STREAM TUBES ( NEWTONS /UN I T H) 

* //( 7( IH 12, E 14.6 ) ) > 

11C FORMAT! 31HOTOTAL THRUST (NEWTONS/UNIT H) E14.6/27H TOTAL POWE 

*R (WATT/UNIT HI E14.6) 

111 F0RMAT!46HCINI TIAL CURRENT DENSITIES (AMPS/lUNIT H)**2) // 

* ( 7! 1H 12, E 14.6) ) ) 

112 FORMAT 1 23HCTOTAL INITIAL CURRENT 3 E12.6, 5H AMPS ) 

113 FORMAT! 3 7H0 TR AN SMI T TED CURRENT AT ACCEL. GRID 3 E12.6.15H AMPS W 
*H ICH IS F6.2, 33H PERCENT OF THE INITIAL CURRENT. ) 

114 FORMAT (37H0TRANSMITTED CURRENT AT DECEL. GRID 3 E12.6.15H AMPS W 
*HICH IS F6.2, 33H PERCENT OF THE INITIAL CURRENT. ) 

115 FORMAT! 25H0INITIAL CURRENTS (AMPS! //(7(1H I2,E14.6))> 

116 FORMAT ( 47H0THRUST DISTRIBUTION BY STREAM TUBES (NEWTONS) // 

* C 7 ( 1 F 12,614.6)1) 

END 
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C TRAJ CALCULATES THE TRAJECTORY CCORUINATES AND VELOCITIES 
SUBROUTINE TRA J ( M ,KE , KEO I 

COMMON AKEM, AT X ( 5 I ) , AT Y I 5 1 ) * AX, AYI51), BASE, CUI51), CUDI51), 

1 DC, DCC , DELY, DX, EPS, ERI40), E T X I 51 ) , ET Y ( 5 I ) , H, I AS , J AS ( 20 ) , JD 

2 , JOT , JTI3C00), KAB, KABB, KAN, KAT , KATT, KATL, KAT2 , KBI13), 

3 KBA, KBB, KBF , KCHI5I), KRL » KTU850), LAST, LB12), LCI3), MO, 

4 NAJ, NPII, NRD, NRL • NT, NT J , NURL, NXEP, PTXI99), PTYI99), 

5 RCUI51), RHI3000), RHDOWN , RHUP, RIN, RX, SAUI5I), SEM, 

6 SIZE, U I 3000 ) , UB ( 3000) , IRH (3000 ) , VA , VAT, VBT , VXI51), 

7 VY (51), XO, XEP , XK, XMP, XQM , XR, XT(1850),YEP 

JE = KE 

JED=KED 

K=M 

AC=A Y ( K ) /D X 
JX = A0 
XA = JX 
XA=AD-XA 
JP= JX+JE 
J S = 4 

JU= JP— JD 

UL=( I .0— XA ) *U( JQ ) + XA*U ( JQ + 1 ) 

UK= i l.C-XA)*U( JP)+XA*U( JP«-1 ) 

C CALCULATION OF LEFTHAND DERIVATIVE 
IF ( XA— . 5 ) 1,1,5 

1 I F < J X ) 2,2,3 

2 YLA=2.0*XA* ( U( JQ + 1 ) — U t JQ) I 
GO TO 4 

3 YL A= ( XA + • 5 ) *U( JQ+ l )-2 . 0* XA* U ( J C ) + (XA-.5 ) *U ( JQ-1 ) 

4 DY=V Y ( K ) / VX ( K ) 

JOX=JX 

GO TO 7 

5 JX = JX + 1 
JQ= JQ+ 1 

X A= XA- 1 .0 
I F ( J X— JO + 1 ) 3,6,6 

6 YLA=2.C*XA*( Ul JQ-1)-U( JO) ) 

GO TO 4 

C CALCULATION OF RIGHTHAND DERIVATIVE 

7 XB = XA +DY 
IF(DY) 9,9,10 

8 JX = J X— 1 
XB= I.O+XB 

9 IF1XB+.54 8,12,12 

11 XB=XB— 1.0 
JX= JX+1 

10 1FIX6-.5) 12,12,11 

12 IF(JX) 14,17,16 

13 JX=JX+2*(1-JD) 

14 JX=- JX 
XB=- XB 

15 JP=JE+JX 

YRA=1 XB+.5 )* U( JP + 1 ) — 2. 0*XB*U(JP) + (XB-.5) *U( JP— 1 ) 

IF(XB) 19,20,20 

16 IFlJX-JD+1) 15,18,13 

17 JP= JE 
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YRA=2.0*XB*lUl JP+1)— U( JP) ) 

XB=A6S( X8 ) 

GO TO 2C 

18 JP=J£+JX 

YR A=2 .C*XB* ( UIJP-1 ) — U t JPJ } 

19 JP=JP-1 
XB=1.C-ABS( XB) 

20 JQ= JP— JO 

USN= ( 1.0-XB) *UI JQ)+XB*U< JG+1) 

UQN= ( l.G-X6)*U< JP1+XB*U( JP+1) 

DUX=. 5* ( UK— UL— USN+UQN) 

VXB = VX ( K )#*2— 2.*XQM*DUX 
VXX=SQRT ( ABS I VXB ) ) 

DT = 2.C*DX/( V XX + VX { K ) ) 

JS=JS-1 

C YA=Y ACCELERATION 

YA=XD*< YLA+YRA ) 

C CY=DELTA Y INCREMENT 

DY=DT*(VY(K)-.5*YA*DT) /OX 
JX=JOX 

IF(JS) 2 1* 21* 7 

21 IFIVXB) 28,29,29 

28 WRITE! 6, 100 )K , JE , JED , JP , JO ,A Y ( K ) , DUX , VX { K ) , V Y { K ) ,VXB,AX 
29 VX1K)=VX* 

VYI K J = V Y (K)-YA*DT 
AY(K)=AY(K1+DY*DX 

C REFLECTION OF TRAJECTORIES IF OUTSIDE BOUNDS 
IF!AY<K))22,23,23 

22 AY ( K ) =-A Y ( K ) 

IF C IAS.LE.O) GO TO 25 
AY ( K ) =— 1 . 

CUI K ) =0 . 

CU (K+l J =0 • 

GO TO 26 

23 IF ( AYI K ) — R I N ) 26,26,24 

24 AY< K )=R IN+RIN-AYIK) 

25 VY ( K ) =- VYI K ) 

27 KC H ( K ) = KCH ( K } + 1 

26 CONTINUE 
RETURN 

100 FORMAT 127H0TRAJ ATTEMPTS TQ TURN BACK / / ( 5 1 5 , 6E 1 5 . 6 ) ) 
END 
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C CORRCT CONDITIONALLY TERMINATES OR SHIFTS THE TRAJECTORIES 
SUBROUTINE CCRRCT (JR) 

COMMON AREM, ATX( 51 ) , ATY( 51 ) , AX, AY151), BASE, CUI51), CUDI51), 

1 DC, DCC, DELY, OX, EPS, ERI4C), E TX ( 51 ) , ET Y I 5 1 1 , H , I AS , JAS I 20 ) , J D 

2 , JOT , JTI3C00), KAB, KABB, KAN, KAT, KATT, K AT 1 , KAT2 , KBI13), 

3 KBA, KBB, KBF , KCHI51), KRL , KTI1850), LAST, LB(2>, LC ( 3 ) , MO, 

4 NAJ, NPIT, NRD , NRL , NT, NT J , NURL, NXEP, PTXI99), PTY199), 

5 RCU (51), RH( 3000) , RHDOWN, RHUP, RIN, KX , SAU(5i), SEM, 

6 SIZE, UI3C00), UBI3000), LRH 1 3000 1 , VA , VAT, VBT , VX(5l), 

7 V Y ( 51) , XO, XEP , XK, XMP, XQM , XR, XT(1850),YEP 
C JN=MESH COLUMN NUMBER 

JN= JR 

IF( ( (JN-KATI’MJN-KATT) J.GT.O) GO TO 63 
J = JN— KAT+1 
GO TO 82 

63 1F( ( ( JN-KAT1)*! JN-KAT2) J.GT.O) GO TO. 87 

J=JN-KATl+KATT-KAT+2 
82 DO 64 K = 1 * NT J 

IFI KCPIK I.LT .0 ) GO TO 64 
C BEGIN IMPINGEMENT CHECK 

IFIAYIK).EQ.-1..0R.AYIK).GT.ER(J) ) GO TO 64 
AA = A Y ( K ) 

IF(K.EQ.l) CU( l) =0. 

IFI K.EQ. l.OR .( K.GT. l.ANO. AY (K-l ) . EQ.-l. ))G0 TO 10 
C CORRECT CU(K) WHEN A Y ( K- 1 ) HAS NOT TERMINATED 
AYD=AY( K-l )-AA 
AYDS = A Y ( K— 1 ) -ER I J ) 

CUI K )=AYDS/AYD*CUC K) 

AYDL=AYO-AYDS 

VYIK ) = < VYI K)*AYDS + VY(K-1)*AY0L) /A YD 
VXIK ) = < VXIK J*AYDS + VX(K-1)*AYDL) / A YD 
AY I K ) = ER ( J ) 

IFIK.GE.NTJ) GO TO 13 

IFIAYIK— II . LE . ER ( J ) .AND.AYI K+ 1 ) .LE.ER(J) ) GC TO 11 
IFI AYIK+1) -GT.ERI J I I GO TO 12 
CUI K+ 1 ) = C. 

GO TO 14 

C CORRECT CUlK+i) WHEN AYIK+1) HAS NOT TERMINATED 

12 AYD=AY I K + l J— AA 
AYDS=AYIK+i)-ERIJ) 

CUI K + l )=AYDS/AYD*CUI K+i ) 

AYDL=AYD-AYDS 
AY I K )=ER I J I 

VY( K ) =IAYDS*VY I K) +AYDL*VY( K+l) I /A YD 
VXIK)=(AYDS*VX(K)+AYDL*VX(K+1) I /A YD 
GO TO 14 

10 IFI K.LT .NTJ .AND.AYIK+l ) .LE.ERI J I ) GO TO 11 
IFIK.EU.NTJ) GO TO 13 

GO TO 12 

11 A Y I K ) =— 1 . 

CU(K+1J=0. 

GO TO 64 

13 IF l ER I J I .GE .RIN) GO TO 11 
C CORRECT FOR LAST TUBE 

AYD = R IN-AA 
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AYDS=R IN-EK ( J ) 

AYDi. = AYC-AYDS 
AY I K ) =ER ( J ) 

CUCK+l )=AYDS/AYD*CUIK+1 ) 
VYIK)=AYDS/A YD*VYIK) 

N=A X/DX •» 1.5 
N=N*LC ( 3 ) 

NN=LC ( 3) 

DIF=ABSI U(NN)-UIN) ) 
V=StlRT(2.*XQM*DIF J 
VX(K)=(AYDS#VX(K)+AYDL*V)/AYD 
14 IF(ER( J ) .GE.RIN) GO TO 11 

64 CONTINUE 

87 RETURN 

END 
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SUBROUTINE PEQ 

COMMON AREM, AT X { 5 11 , ATY ( 5 I ) , AX, AYI51), BASE, CU(51), CUD(51>, 

1 DC, CCC, DELY, DX, EPS, ER(4C), E TX ( 5 1 ) , ET Y ( 51 ) , H , I AS , JAS < 20 ) , J D 

2 , JOT , JT13000), KAB, KABB, KAN, KAT , KATT, KAT 1 , KAT2 , KB(13), 

3 KBA, KBB, KBF , KCH(51), KRL, KTU850), LAST, LB<2>, LC(3), MO, 

4 NAJ, NP IT, NRD , NRL , NT, NT J , NURL, NXEP, PTX{99), PTY(99), 

5 RCU (51), RH ( 3C00 ) , RHDO InN , RHUP, RIN, RX , SAU(51J, SEM, 

6 SIZE, U(3C00i, UB { 30001 , LRH (3000) , VA , VAT, VBT , VX(51), 

7 V Y ( 51) , XD, XEP , XK, XMP, XQM , XR, XT(1850),YEP 

C PEQ CALCULATES THE EQUIPOTENT IAL LINE FOR THE CU CALCULATION 
PO T EN = U ( KBA ) 

DX=H 

AL= . 1*H 
NO=KABB+ 1 
AX = 0 .0 
L = 0 

JE=LC( 2) 

JD=LC 13) 

20 DO 8 J J= 1 , KAB 
AAY = 0 . 0 
JED = JE -*JD- I 
DO 7 K= J E, JED 
J=K- JD 

I F ( (U<KI-POTEN)*CU< J)-POTEN) ) 1,1,7 

1 DIF=ABS(U(J)-U(K) ) 

L=L + 1 

C COORDINATES ARE TAKEN FOR EQUAL DELTA Y ONLY 
PTYI L )=AAY 
IF(DIF) 2,6,2 

2 DO 15 KJ=1,NG 

IF( ABSI AAY-ATY ( KJ ) )-AL) 16,16,15 

15 CONTINUE 
GO TO 4 

16 IFl AX-ATXIKJ ) ) 3,3,4 

3 HK=AX+DX— ATX(KJ ) 

AK = A TX ( K J ) 

GO TO 5 

4 HK = DX 
AK = AX 

5 PTX IL )=ABS ( U( J ) -POTEN) /DIF*HK+AK 
GO TO 7 

6 PTX ( L ) =A X+DX 

7 AA Y =AA Y +DX 
AX=AX +DX 

8 JE=JE+JD 

C SORT Y-COORD IN INCREASING Y-ORDER 

DO 1 1 J= 1, L 
LL =1 -J + l 
T=0 .0 

DO 10 1 = 1, LL 
IFCT— PTYIIJI 9,9,10 

9 T = P T Y I I ) 

NN= I 

10 CONTINUE 
PP=PTY ( LL ) 
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PT Y { LL J - T 

P T Y ( NN )=PP 
PP=PTX(LU 
P TX ( L L ) =P TX ( NN ) 

PTXINN) =PP 

11 CONTINUE 
KAP^O 

C C I SCAR D DUPLICATES 

J-2 

31 1FIPTYIJ )— PTY(J-l) ) 14,12,14 

12 NN=L-1-KAP 
KAP=KAP+1 

DO 13 J J = J , NN 
PTXI JJ ) = PT X { J J+ 1 ) 

13 PTYI JJ J=PTYIJJ+1) 

14 J=J+1 

1FI J. Lb. ( L-KAP ) ) GO TO 31 
IFIJOT.LE.O) GO TO 30 
WRITE (6, 100 ) PO TEN 
NAL = L-KAP 

WR ITE <6, 101 ) ( J ,PTX ( J ) , PT V ( J ) , J=1 , NAL ) 

30 RETURN 

100 FORMAT { 57HOCURRENT DENSITIES ARE CALCULATED USING E QU I P OT ENT I AL 

* OF F 1 0 * 5 , 33H VOLTS WHICH HAS X-Y COORDINATES ) 

101 FORMAT ( 7 ( 1H 12, 2H ( F 5 . 3 , 1 H , F5 . 3 , 2H ) )) 

END 
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INTO EQUAL ARC LENGTHS 


C ARC DIVIDES TEE EMITTER 
SUBROUTINE ARC 

COMMON AKEM , A TX ( 5 I ) , AT V < 5 1 ) , AX, AY(51), BASE, CU(51>, CUD(5I), 

1 DC, DCC , DELY, DX, EPS, ERC4G), E TX ( 5 1 > , ET Y( 5 1) , H , I AS , J A S ( 20 ) , J D 

2 , JOT , JTC3C00J, KAB, KABB, KAN, KAT, KATT, K AT I , KAT2 , KB(13), 

3 KBA, KBB , KBF , KCHI51), KRL , KTC1850), LAST, LB(2), LC { 3 ) , MO, 

4 NAJ, NPIT, NRD , NRL , NT, NT J , NURL , NXEP, PTX(99), PTY<99), 

5 RCU< 51), RH( 30CC ) , RHDOWN , RHUP, RIN, RX, S AU ( 5 1 ) » SEM, 

6 SIZE, UI3C00), UBC3000), LRH ( 3 000 ) , V A , VAT, VBT , VX(51), 

7 V Y ( 5 1 J , XC, XEP, XK, XMP, XQM, XR, XT(1850),YEP 

DIMENSION LQ(ICO) 

EQUIVALENCE! LQ( I) ,RH( 1) > 

C K ABB = NUMB ER OF INPUT COORDINATES FOR EMITTER LESS ONE 
NBH=N T J + I 
B h=NBH 
SEM = 0 .0 

C SUMMING TOTAL ARC LENGTH 
DO I J- 1 , KABB 

SA= SQR T { (ATX { J+l )-ATX( J) ) **2 + { AT Y ( J + 1 )-ATY( J) )**2 ) 

1 SEM=S£M+SA 

C AR EM = DELT A ARC LENGTH 

AR EM = SEM /BH 
K= I 

R EM = 0 .0 
J = I 

HYP =AREM 

C CALCULATION OF COORDINATES 
ETX! K)=ATX(K) 

ETY I K ) = A TY ! K ) 

K = K + 1 

2 XDIF=ATX( J )-ATX( J+I) 

YDIF=ATY(J+1 )-ATY! J) 

SA=SQRT( XDIF**2+YDIF**2) 

C S A = YD IF/SA 

SNA = XQ IF/SA 

3 ETX(K)=ATX!J )-HYP*SNA 
ETY ( K ) = A T Y < J ) +H YP *C S A 
HYP = AR EM +HYP 

K=K + 1 

R EM = R EM +AREM 
IF(K-NBH) 4.4,6 

4 IF! SA-REM-AREM) 5,5,3 

5 J=J+1 
HYP=AREM-SA+REM 
REM=HYP- AREM 

GO TO 2 

6 CONTINUE 

ETXCK ) = ATX(KABB+i ) 

ETYI K3=ATY(KABB+I) 

C CHECK ON LOWER SYMMETRY 
L2 IF ( LAST ) 17,18,18 

17 K = K + 1 

NT J = N T J + 1 

XDIF = ETXIK- D-ETX (K-2) 

YOIF = ETY<K- 1 HETY! K-2) 
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SA=ATAN(XDIF/YDIF ) 

ETX { K )=ETX( K-l ) +AKEM*SI N( SA) 

ETY (K)=ETY(K-1)+AREM*C0S( SA) 

C CONDITIONAL PRINTOUT 
18 IF(JOT) 8,8,7 

7 WRITE ( 6 , 103 )SEM,AREM 
NAL=KABB +1 

WRITE (6*100 )(J,ATX(J) , A T Y ( J } ,J = 1 ,NAL ) 

NAL =K- 1 

IF(KBF) 13, 14, 14 

14 K = 1 

GO TO 15 
13 K = 2 

15 NUM = 1 

CO 16 J =K , ICC 
LOU ) =NUM 

16 NUM =N LM + 1 

WRITE (6,102) ( LQ ( J ) , E T X ( J ) ,E T Y (J ) , J= K , NA L ) 

DO 20 J=1,10C 
20 RH( J ) =0 • 

8 RETURN 

ICO FURMAT { 24H0X- Y EMITTER COORDINATES // (7< 1H 12, 2H ( F 5. 3, i H, F5 . 3, 

* 2H ) ) } 1 

102 FORMAT ( 2 8H0X-Y BEGIN TRAJ. COORDINATES / / ( 7 C 1 H 12, 2H (F5.3, 

* 1H,F5.3,2H> ))) 

1 C 3 FORMAT (4 1H0EMI TTER ARC LENGTH , DELTA ARC LENGTH 2F10.5) 

END 
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C RTEST CHECKS ON THE UPPER AND LOWER BOUNDS 
SUBROUTINE RTEST 

COMMON AREM, ATX { 51 ) ,A TY ( 5 1 ) , AX, AY<51), BASE, CU<51), CUD(51), 

L DC, DCC , DELY, DX, EPS, ER140), E TX ( 5 I } , ET Y { 5 l ) , H , I AS , J AS ( 20 I , JD 

2 , JOT , JTI3C00), KAB, KABB, KAN, KAT, KATT, KAT1 , KAT2 , KB(13), 

3 KBA, KBB » KBF , KCHI51), KRL , KTU.850), LAST, L8(2), LC<3), MO, 

4 NAJ, NPIT, NRD , NRL , NT, NT J , NURL, NXEP, PTXI99), PTY(99>, 

5 KCUI51), RHI3C00), RHDQWN , RHUP, RIN, RX, SAU<51), SEM, 

6 SIZE, U I 3000 ) , UBI3000), CRH 1 3000 ), VA , VAT, VBT , VX(5I>, 

7 VY151), XD, XEP , XK, XMP, XQM , XR, XT(1850),YEP 
C KAN=TEST POINT 

IFI KRL.EQ.NRL) WRITEI6.102) 

RT=RH ( KAN ) 

SX=1.C 

C NO CHECK IF CURRENT DENSITY IS SPECIFIED 
IFIKBB )25, 24,25 

25 IF( KRL ) 19, 26 » 1 9 

26 IF ( MO i 11,19,11 

24 IF ( KRL ) 22, 23,22 

23 IF ( MO ) 11,22,11 

C FIRST TWO CYCLES SET UPPER AND LOWER BCUNDS 
22 IF ( KRL-NRL+ 1 ) 5,2,1 

C RHUP=UPPER BOUND 

1 RHUP=RT 
GG TO 19 

2 IF ( RT— RHUP ) 18,18,3 
C RHDOWN = LOW ER BOUND 

3 RHDOWN=RHUP 
MO=-MO 
RHUP=RT 

4 GO TO 19 

C M0=+ 1 .CHECK IS ON UPPER BOUND 
C MQ=-1 .CHECK IS ON LOWER BOUND 
C MO =0 , NO CHECK 

5 IF1MC) 10,19, 8 

6 SX=SX#RX 

DO 7 J = 1 , NT 

7 RH( J )=RH( J )*RX 

8 IF(RFUP-RHIKAN) ) 6,9,9 

9 IF{RH(KAN)-RHDOWN ) 11,16,16 

10 I F ( ( RT-RHDOWN ) << (RT-RHUP) ) 16,16,11 

11 P0W=1.0/SX 
CO 12 J= 1 ,NT 

12 RH( J )=.5*(RH( J I^POW+URH ( J) ) 

KRL = 1 

J J = NKL * 1 
KB ( J J ) = 7 7 
KB( JJ+1 )=77 
WRITE 16,100) 

IFIKBB.GT.O) GO TO 15 
IF ( MO ) 14,13,13 

13 RHUP =R T 
GO TO 15 

14 RHDOwN=RT 

15 M0=0 
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16 IF { MG } 18,19,17 

17 RHUP =Rh 1 KAN J 
GO TO 19 

18 RHDGUN = R T 

19 XT { 1 ) = .25*R'HOOWN 
XT 1 2 ) = * 2 5*RH f K AN ) 

XT { 3}=.25*RHUP 

WR 1 T E ( 6 t 10 1 ) X T ( 1 ) , K AN » X T ( 2 ) , XT(3) »KAN ,U ( KAN ) 

I NRL-KRL + 1 
IF(KB(IV»)J 21,21,20 
C TROUT PRINTS OUT THE RHS 

20 CALL TROUT 

21 MO=-MG 
KRL=KRL- 1 
RETURN 

100 FURMAT1 11H0RP=AVERAGE ) 

101 FORMAT ( 8HORHDOWN=F1 1 • 7 ,4H RH ( I 4 , 2H ) =F 1 1 . 7 ,6H RHUP=F1 

* 2HJ=F15.7J 

102 FORMAT ( 26H0START CF POISSON SOLUTIGN ) 

END 


• 7 , 3 H U< 14, 
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C TRCUT PRINTS CL T THE RHS 
SUBROUTINE TROUT 

COMMON AREM, ATX ! 5 1 ) , AT Y < 5 1 > , AX, AY(51), BASE, CU!5i), CUD(51), 

1 OC, DCC, DELY, DX, EPS, ERIAG), E TX ( 5 1 ) , ET Y ( 5 1 I , H , I AS , J AS ( 20 > , J 0 

2 , JOT , JTI3C00), KAB, KABB, KAN, KAT , KATT , KAT1, KAT2 , KB(13), 

3 KBA, KB8, KBF , KCH151), KRL » KTI1850), LAST, LB12), LC ( 3 ) , MO, 

A NAJ, NP IT » NRO , NRL, NT, NT J , NURL, NXEP, PTX!99), PTY<99), 

5 RCU (51), RH( 3000 ) , RHDOWN , RHUP , RIN, RX, SAUJ51), SEM, 

6 SIZE, UI3C00), UB (3000) , LRH ( 3000 J , VA , VAT, VBT , VX151), 

7 VY 151), XD, XEP, XK, XMP, XQM » XR, XT( 1850),YEP 

IF ( MO ) 1,8,1 

C NORMAL PRINT OUT 

1 WRITE (6, 100) NT 

WRI TE (6, 10 A ) 

J=1 

DO 3 K = 1 * NT 

XT ( J )= .2 5*RH I K ) 

KT ( J ) = K 
J=J + 1 

I F { J — 9 ) 3,2,2 

2 WRITE !6, 101 ) ( KT(M) ,M = 1 ,8) , (XT! M) ,M = 1 ,8 ) 

J = 1 

3 CONTINUE 
IFIJ-2) 7, A, A 

A 00 5 K = J * 8 

KT ( K )=0 

5 XT ( K ) = 0. 0 

WRITE I6*1Q1)(KT!M)»M=1,8) , ( XT { M ) ,M = 1 ,8) 

7 RETURN 

C FINAL PRINT CUT IS AVERAGE 

8 WRITE 16,102) 

WRITEI6, 1031 

17 WRITE 16, ICO) NT 
WR I T E ( 6 , 10 A ) 

J = 1 

CO 10 K= 1 , NT 
K T ( J ) =K 

XT l J ) = . 125*1 RHiKH-URHIK) ) 

J=J + 1 

IF(J-S) 10,9,5 

9 WRITE ( 6,101 ) ( KT(M) ,M = 1 ,8) , ( XT! M) ,M=1 ,8) 

J=1 

10 CONTINUE 
IFIJ-2 ) 7,11,11 

11 DO 12 K= J , 8 
KT I K ) =C 

1 2 XT i K ) = C • C 

WRITE 16,101 IIKT(M) ,M = 1,8) , IXTIM) ,M=1 ,8) 

100 FORMAT! 1FCI5, 10F RH VALUES) 

101 FORMAT! 1H 8I5.8F11.A) 

102 FORMAT (A7H0RHQLTPLT IS AVERAGE OF THIS AND PREVIOUS CYCLE ) 
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103 FORMAT ( 27H0C ON VERGED POISSON SOLUTION ) 

104 FORMAT ( 1FG, 1CX, 2 OH MESH POINT NUMBERS ,49X,10H RH VALUES 
ENO 
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C TRCU CALCULATES THE CURRENT AT THE EMITTER AND INITIALIZES 
C THE TRAJECTORIES 

SUBROUTINE TRCU(KE,KED) 

COMMON AREM, AT X { 5 1 ) , A TV { 5 1 ) , AX, AYI51), BASE, CU{51), CUDI51), 

1 CC, CCC, DELY, DXt EPS, ER(4C), E TX ( 51 ) , ET Y t 5 1 ) , H , I AS , JAS i 20 ) , J D 

2 , JOT , J T ( 3C00 ) , KAB, KABB , KAN, KAT, KATT , K AT I , KAT2 , K B ( i 3 ) , 

3 KBA, KBB , KBF , KCHI51), KRL, KTC1850), LAST, LB(2), LC { 3 ) , NO, 

4 NAJ, NP IT, NRD , NRL , NT, NT J , NURL, NXEP, PTX<99>, PTY(99), 

5 RCU (51), RH( 3000) , RHDGWN, RHUP , RIN, RX, SAUCSUt SEM, 

6 SIZE, U { 3 COO ) , UBI3000), LRH 1 3000 ) , V A , VAT, VBT , VX(5L), 

7 V Y ( 5 1 ) , XC, XEP, XK, XMP, XOM, XR, XT(1850),Y£P 
DIMENSION U S ( 2) ,UR (2) 

XU=ABSI XOM) 

JE = KE 
J ED = KED 
AL A = • 2*H 
AC A N =1.0 
ISW = 0 
K= I 

116 IS=C 

LL = 0 
KK = K 

C CHECK ON TERMINATED TRAJECTORIES 
IF(AYCK)+1.0) 16,16,18 
C CHfcCK ON UN IN I AL I ZE D TRAJECTORIES 

18 IF(KCHIK)) 20,15,15 

20 IF I 1 Sw ) 22,22, 1 

22 I F ( AY(K)*ACAN) 1,1,19 

19 I S= 1 

I Sto = 1 

C CHECK IF TRAJECTORY CAN BE INIALIZED 

1 IFIAX-ETX(K+1) -ALA ) 21,21,2 

2 I F ( AX— ETX(K+2 ) - ALA ) 21,21,5 

5 N=— 1 

AC AN = 1 . C 
KK = KK + 1 

TA= i ETX ( KK- 1 J-ETX IKK) ) / IETV IKK-l )-E TY IKK) ) 

TB= ( E T X { KK )-ETX{ KK + 1) ) / ( ET Y { KK ) -E T Y ( K K + 1 ) ) 

Tb=-.5*(TA+TB) 

DE = A TAM TB ) 

33 XX=ETY|KK)-TB*ETX(KK) 

J = 1 

3 N = N + 1 

4 J J = J + N 

TA=PTY{ JJ)-PTY(JJ+1) 

TC=PTX( JJ)-PTX{ JJ+1) 

YY=TC*PTY( JJ )— T A* P T X ( J J ) 

CPXM TC*XX-YY) / { TA-T6*TC) 

CPY=XX+CPX*TB 

7 IF (PTYCJJ+l )-CPY) 3,8,8 

8 LL = LL + 1 

24 I F ( AY(KK-l) )9, 9, 10 

C TRAJECTORIES INITIALIZED 

9 AY I KK- 1 ) =E TY(KK)+TB*(AX-ETX (KK) > 

ACA N = 0 .0 
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XA = A Y { KK— 1 ) / UX 
JX = X A 
JP = J £ + J X 
AD = J X 
XA=XA-AD 

UP = ( l.C-XA)*U( JP)+XA*Ui JP+l) 

V = SQRT< 2.0*XQ* (ABS( VA~UP) ) ) 
VXIKK-1)=V*C0S{DE) 
VYIKK-1}=V*SIN(DE) 

10 US l LL I =CPX 
UR( LL }=CPY 
IF(K-NTJ) 12,11,12 

C C SPECIAL FOR LAST CURRENT DENSITY 

11 HT = AR EM +AREM 

IFCLASTJ 40,41,41 

40 HT = 0 . C 
GO TO 43 

41 KK=KK + 1 
LL= 2 

29 US 121 = PTXCJJ+l) 

UR 1 2 ) = PTYIJJ+1) 

PPX=.5* (ETXJ KK~1 J+ETXCKK) ) 

P X = • 5# I US f 1 ) +US ( 2 ) ) 

PPY = .5* (ETYIKK-1 KETY(KK) ) 

P Y= • 5* ( UR C 1 ) +UR ( 2 ) ) 

GO TO 31 

12 IF ( LL“ 2 ) 5,13,13 

13 HT =AR EM 

23 PPX=.5*(ETX(KK~1J+ETX{KK) ) 

PX= • 5*( U S ( 1J+US(2 ) ) 

PP Y = . 5* { ETY { KK-1I+ETYIKK) ) 
PY=.5*(URll)+URI 2) ) 

31 DEX = SORT { { PX-PPX) **2+(PY-PPY)**2 ) 

C CALCULATION OF CURRENT DENSITY 

17 CU{ K+l )=HT/DEX**2*XK 

CUDIK+1)=CU(K+1)/HT 
I F ( KBB 42 , 42 , 4 3 

43 CUIK + 1 )=RCU(K+1I*HT 

CUDIK + 1 ) =RC U l K+ 1 ) 

42 SAU(K*1)=CUIK+1) 

KC H { K ) =0 

I F ( K- 1 ) 14,14,26 

C SPECIAL FOR FIRST CURRENT DENSITY 

14 D£ = E TX ( 1 ) — E T X I 2 ) 

102 DE=-CE/(ETY ( D — ETY12) ) 

104 XX =£ T Y ( 1 )-DE* ETXI 1) 

J =0 

110 J = J + 1 

109 TA=PTY (J )-PTY( J+l) 

TC = P TX ( J )-PT X ( J + 1 ) 

YY=TC*PTY( J )~TA*PTX< J) 

CPX=( TC*XX-YYJ /( TA-TC*DEI 
CP Y=XX+CPX*DE 

107 IF(PTY< J + D-CPY) 110,108,108 

10 8 PX=.5*<ETXm+ETXI 2) ) 

PPX=.5*(CPX+US( 1 ) ) 
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Ill 

113 


114 


45 

44 


118 

115 

47 

46 

26 

21 


PY=.5*i E TY( 1 ) *E T Y { 2 ) ) 
PPY=.5*(CPY+UR( 1) ) 

DEX = SQRT t {PX-PPX>**2+( PY-PPY) **2) 

HT = 2 . 0*AREM 

IF(KBF) 115,115,113 

HT=.5*HT 

NTJ=NTJ+1 

NAJ=NAJ4 1 

J R = N T J 

DO 1 1 4KL =2 , NT J 
SAU ( JR+l)=SAUt JR) 

CUD t JR + 1 ) = C UD ( JR ) 

RCU<JR+1 )=RCU t JR) 

CU( JR+l)=CUt JR) 

AY t JR ) =A Y ( JR- 1 ) 

VY( JR ) = V Y { J R- 1 ) 

VX< JR)=VX( JR-1) 

KCH t JR ) =KC H ( JR-1) 

JR = JR— 1 

RCU ( 2 ) =RCU ( 1 ) 

R CU t 1 ) =G • 0 
CUD t 1 )=0 .0 
CU( l)=C.O 

cut 2)=HT*XK/DEX**2 
CUD (2 )=CU( 2)/HT 
I F { KBB ) 4 4 , 44 ,45 
CU{ 2 )=HT*RCU!2) 

CUD ( 2 ) =RCU ( 2) 

SAU l l)=CUt 1) 

SAU < 2 ) =C Ut 2 ) 

A Y ( 1 ) = E T Y ( 1 )+DE* t AX-ETX( 1 ) ) 
TB=ATAN1DE ) 

XA = A Y ( 1)/DX 
J X = X A 
JP= JE+JX 
AD= J X 
XA = XA— A D 

UP = ( 1.0-XA)*Ut JP) + XA*Ul JP+i) 
V=SQRT(2.C*XU*A6S( VA-UP) ) 

V X 1 1)=V*C0SITB) 

V Y ( 1)=V*SIN(TB) 

K = K + 1 
J= 51 

DO 118 1=1,50 
ETXt J) = ETX( J-l) 

ETY ( J ) = ETY t J- 1 ) 

J=J-1 
GO TO 26 

CUt 1 )=HT*XK/DEX**2 
CUD t 1KU 1)/HT 
IFt KBB )46, 46,47 
cut 1 )=HT*RCUt 1) 

CUD t 1 ) = RCU ( 1 ) 

SAUt 1)=CU( 1 ) 

CONTINUE 

IFt IS) 16,16,15 


C TRAJ CALCULATES THE TRAJECTORIES AFTER THEY HAVE BEEN INIALIZED 

15 CALL TRAJ(K, JE, JED) 

16 K=K + 1 

IF(NTJ-K) 117,116,116 
117 RETURN 

END 
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SUBROUTINE C AL R < KE f KED ) 

COMMON AREM, A TX < 5 1 ) , A T Y ( 5 1 ) , AX, AY(51), BASE, CU(51), CUD ( 5 1 ) , 

1 CC, DCC , DELY, DX, EPS, ER(4C), E TX ( 51) , E T Y { 5 1 ) , H , I AS , J AS ( 20 ) , JD 

2 , JOT , J T ( 3 CUO ) , KAB, KABB, KAN, KAT , KATT, KAT I , KAT 2 , KB(I3), 

3 KBA, KBB • K6F , KCH(51), KRL , KTU850), LAST, LB(2), LC(3), MO, 

A NAJ, NPIT, NRD, NRL , NT, NT J , NURL, NXEP, PTX(99), PTY(99), 

5 RCUl 51), RH ( 300C ) , RHDOWN, RHUP, RIN, RX, SAUI5L), SEM, 

6 SIZE, U { 3C00 ) , UB ( 3000 ) , LRH < 3 000 ) , VA , VAT, VBT , VX(5L), 

7 V Y ( 5 1 ) , XD, XEP , XK, XMP, XQM, XR, XT(1850),YEP 

JE=KE 
J EU = KED 
SWA=0.0 

DO 33 J J = 1 , NAJ 
J0 = 0 
J = J J 

C I F CU ( J ) =0 * , NO CONTRIBUTION 
IF(CUIJ)) 33,33,1 
C IF CU( U, SPECIAL 

1 I F ( J -1) 2,2,5 

2 IF(KChtJ)) A, A, 3 

3 IF(VYU)) 35, A, A 

3 5 HT=2«0*(2.C*RIN-AY{J ) ) 

NA = J£ 

NB= JED 
WA=VX ( J ) 

W Q = W A 
JG= 1 

I F ( SWA ) 39, 39, AC 
39 SWA = 1.0 

XH = AY 1 J )-2.0*RIN 
YL=-XH 
GO TO 23 
AO SWA=O.C 

XH= A Y ( J ) 

YL=XH*HT 
GO TO 23 

A HT=2.0*AY(J) 

XH=- A Y < J ) 

J X = A Y ( J ) /DEL Y+ 1 .0 
NA = JE 
NB= JX+JE 
WA = VX ( J ) 

WB= W A 
YL = A Y ( J ) 

GO TO 23 

C IF CU( NAJ ), SPEC IAL 

5 IF(J-NAJ) 9,6,33 

6 IF(KCMJ-l)) 8,8,7 

7 J = J-1 

IFIVYC JJ-1)) 8,22,22 

8 FT = 2. 0*1 KIN-AY (JJ-1) ) 

XH=AY( JJ-1) 

YL=XH+HT 

JX=AY( JJ-1J/DELY 
NA = J X + JE 
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NB = JED 
WA=VX( JJ-1) 

WB= WA 
GO TO 23 

9 IF(KChtJ)) 11,34,34 

C IF TRAJECTORIES HAVE REF LEC TE D , ADD IN MIRROR IMAGE 
34 IF(KCH( J-1)-KCH( J) > 10,11,15 

10 J0=1 

lF(VYIJ)) 19,16,16 

11 HT=AY<J)-AY( J-l) 

IF(HT) 12,12,13 

12 HT=-H T 

XH= A Y ( J } 

JX=XH/DELY 

NB=AY( J-D/DELY+1.0 
WA=VX { J ) 

wb=vx< j-n 

GO TO 14 

13 XH = A Y ( J- 1 ) 

JX=Xh/DELY 

WA=VX{ J-l) 

WB = VX { J ) 

NB=AY(J)/DELY+1.0 

14 NA = JX + JE 
NB = NB + JE 
YL=XH+HT 
GO TO 23 

15 J0=1 

IFtVY(J-l)) 19,19,16 

16 HT=AY{ JJ+AYl J-l) 

NA = JE 

I Ft SWA ) 17,17,18 

17 XH=— A Y ( J— 1 ) 

YL=AYl J 3 
JX=YL/DELY+1.0 
NB = JX JE 

WA = VX ( J- 1 ) 

WB=VX ( J ) 

SWA = 1 .0 
GO TO 23 

18 SWA =0.0 
XH=— AY ( J ) 

YL=AY( J-l) 

JX=YL/DELY4-1 .0 
NB=JX+JE 

WA = VX I J ) 

WB=VX t J-l) 

GU TO 23 

19 HT=RIN+RIN-AY{ J)-AY< J-l) 

NB = J ED 

I Ft SWA ) 20,20,21 
20 SWA= 1.0 

XH = AY( J ) 

Yi_ =R IN + RIN-AYt J-l) 

JX=XH/DELY 
NA= JX+JE 
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WA=VXt J) 

WB=VXI J-l) 

GO TG 23 
21 SWA = Q . 0 

XH=AY( J-l) 

YL=RIN+RIN-AY{ J) 

JX=XF/D£LY 
NA=JX+JE 
WA=VXI J-l) 

WB=VX{ J) 

GO TO 23 

22 FT=2.0*IRIN+AY(J) ) 

JO= 1 
NA = JE 
NB = J ED 
WA = VX I J } 

WB = WA 

IF1SWA) 3 7,37,36 

37 SWA= 1.0 
XH=- A Y { J ) 

YL=XH*HT 
GO TG 23 

38 SWA =0 . 0 
YL = AY { J ) 

XH= YL— HT 

23 DO 32 K=NA,NB 

XA=K— JE 

IFiXA.LE.O.) JU= J— 1 

25 XUD= ( XA- . 5 ) #DEL Y 

26 XX = XUD +DELY 
YU=AMAX 1 (XUD ,XH ) 

Y0=AMIN1IXX,YL) 

IF 1 YD- YU ) 32,32,27 

27 XA=.5*( YD+YU) 

W = W A + ( XA— XH)*(WB— WA) /HT 
IF ( K- JED ) 29,28,28 

28 JU=NA J- J J 

29 IF(JU) 31,31,30 

30 IF ( JO ) 36,36,31 

36 w=.5*W 

C CALCULATION OF RHS 

31 RH(K)=RH(K)+(YD-YU)*CU(JJ) /(HT*W> 

32 JU=0 

1F1SWA) 23,33,1 

33 CONTINUE 
RETURN 
END 
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ATRCU CALCULATES THE TUBE CURRENTS ANC CALCULATES 

THE TRAJECTORY CUORD AND VELOCITIES IN SUB TRAJ 
SUBROUTINE ATRCU IKE.KED) 

COMMON AREM, ATX ( 5 1 ) , AT Y I 51 > , AX, AYI51), BASE, CUI51), CUDI5I), 

1 DC, DCC , DELY, DX, EPS, ERI4C), E TX t 51 > , ET Y I 51 ) ,H , I AS , JAS ( 20 ) , J D 

2 , JOT , J T < 3000 ) , KAB, KABB , KAN, KAT, KATT, KATl, KAT2 , KBI13), 

3 KBA, KBB , KBF , KCH ( 5 1 ) , KRL , KTI1850), LAST, LB(2), LC ( 3 ) , MO, 

4 NAJ, NPIT, NRD , NRL , NT, NT J , NURL, NXEP, PTXI99), PTYI99), 

5 RCUI51), RH ( 3000 ) , RHOQWN , RHUP, RIN, RX, SAUI51), SEM, 

6 SIZE, UI3000), UB I 3000 ) , LRH ( 3000 ) , VA , VAT, VBT , VXI51), 

7 VYI5U, XD, XEP , XK, XMP, XQM , XR, XT( 1850),YEP 
DIMENS ION US( 2) ,UR 12) 

JE = KE 
JED=KED 
ALA=. 2*H 
ACAN = 1 .0 
XQ=ABSI XQM ) 

ISW = 0 
K= 1 

116 IS=C 

LL = C 
KK = K 

C CHECK ON TERMINATED TRAJECTORIES 
IF ( AY l K ) +1 .0 ) 16. IE, 18 
C CHECK ON UNINIALIZED TRAJECTORIES 


18 

I F ( KCH C K ) ) 20 t 15,15 

20 

IFilSVi) 22,22,1 

22 

I F ( AY(K)*ACAN) 1,1,19 

19 

IS=1 


1 SW=1 


C CHECK IF TRAJECTORY CAN BE INIALIZED 

1 I F t AX-ETX{K + 1)-ALA ) 21,21,2 

2 IFIAX— ETXIK+2)— ALA) 21,21,5 

5 N=— 1 

ACAN = 1.0 
KK = KK+ 1 

TA= ( ET X I KK— l)-ETX(KK) ) / ( ETY ( KK-i )-E TY l KK ) ) 
TB= I ETX(KK)-ETX (KK+1) ) / ( ET Y ( KK ) -E T Y I KK+ 1 ) ) 
TB=-.5*( TA+TB) 

DE = ATANITB ) 

33 XX=ETY(KK)-TB4ETX(KK) 

J=1 

3 N=N + 1 

4 JJ=J+N 

TA=PTYI JJ)-PTY( JJ+1) 

TC=PTX( JJ)-PTX ( JJ+1 ) 
YY=TC*PTYIJJ)-TA*PTX(JJ) 

CPX=I TC*XX-YY ) / ( TA-T8*TC) 

CPY=XX+CPX*TB 

7 IFIPTYI JJ+1J-CPY) 3,8,8 

8 LL = LL+ 1 

24 IF 1 AY I KK— 1 ) >9,9,10 

C TRAJECTORIES INITIALIZED 

9 AY(KK-1>=ETY(KK)+TB*(AX~ETX(KK) ) 

ACAN = 0 .0 
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XA=AY( KK-li/DX 

JX=XA 

AD = J X 

JP= JE+JX 

XA = X A—AD 

UP= ( 1.0-XA)*U( JP >+XA*U( JP+1 > 
V=SQRT<2.C*XQ*ABS< VA-UP) ) 
VX(KK-l>=V*COS<DE ) 

V Y { KK-1)=V*SIN(0E) 

10 US ( LL ) =CPX 
URC LL a=CPY 

IF ( K—NT J ) 12*11,12 

C C SPECIAL FOR LAST CURRENT DENSITY 

11 HT = AR EM +AREM 

IF ( LA ST ) AO. 41,41 

40 RR = 0 *0 
GO TO 45 

41 KK =KK * 1 
LL = 2 

29 US ( 2 ) = P TX ( J J+ 1 ) 

UR ( 2 ) = PTYIJJ+1) 

PPX = • 5* (ETXIKK— 1 )+ETX{KK) ) 

P X = .5<MUS( 1 ) +U S ( 2 ) ) 

PPY=.5*(ET Y(KK-1)+ETY{ KK) ) 

P Y = . 5* ( UR ( 1 ) +UR ( 2 ) ) 

GO TO 31 

12 IF { LL- 2 ) 5,13,13 

13 HT = ARE M 

23 PPX=.5* i ETX( KK— 1 ) +ETX ( KK ) ) 

PX= * 5*( USI l)+US(2)) 

PPY=.5*(ETY{ KK-1)+ETY< KK) ) 

P Y= • 5*( UR { 1 ) +UR i 2 ) ) 

31 UEX = SQRT t (PX-PPX) **2 + i PY-PPY>**2 ) 

C CALCULATION OF CURRENT DENSITY 
17 RR = B A S E-PP Y 

RR=6#2831 8 53*RR 

CU( K + l )= ARE M/DE X**2*XK*RR 

CUD(K*11=CUIK+1 )/( RR*AREM) 

I F ( KBB 2 4 4, 44,45 

45 CU( K + 1 )=RCU( K + 1 ) *RR*AREM 

CUD ( K -t* 1 J =RCU { K + 1 ) 

4 4 SAUIK* W=CUIK+1) 

KCHIK ) =0 

IFIK-1) 14,14,26 

C SPECIAL FOR FIRST CURRENT DENSITY 

14 0£ = ETX( 1 J-ETXI 2) 

102 CE=-CE/IETY< 1)-ETYC2) ) 

104 XX =E T Y ( 1 )-DE*ETXl 1) 

J = 0 

110 J=J+1 

109 TA = PTY( J )-PTY( J+l) 

TC = P TX ( J )— PTX( J + 1 ) 
YY=TC*PTY(JJ-TA*PTXIJ) 

CPX={ TC*XX-YY) /( TA-TC*OE) 
CPY=XX+CPX^DE 

107 IFtPTYt J + U-CPY) 110, ICE, 108 
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108 PX=.5*(ETX( 1) +ETX12) ) 

PPX=.5*1CPX*US ID) 

PY= . 5*( ETY I D+ETYI2D 
RR=BASE~PY 
RR=6.2831853*RR 
PPY = .5*lCPY+liR( 1) ) 

D£X=SQRT ( ( PX-PPX)**2+{PY-PPY) **2 ) 

NTJ=NTJ+1 
NAJ = NAJ+ 1 
JR = N T J 

DO 1 14KL=2 »NT J 
SAUJ JR+l)=SAU( JR) 

RCU ( JR+1 )=RCU( JR ) 

CUD( JR + 1 )=CUD ( JR J 
CU( JR+l)=CUi JR) 

AY( JR j=AY( JR-1 ) 

VY( JR )=VYl JR-1 ) 

VX1 JR)=VX{ JR-1) 

KCH( JR)=KCH( JR-1 ) 

114 JR= JR— 1 

RCU ( 2 ) =RCU I 1 ) 

RCU J 1 )=0.0 
CUD 11 ) = G .0 
CU( 1 )=C .0 

CU( 2 )=AREM*XK/DEX**2*RR 
CUD<2)=CU(2)/(AREM*RR) 

IF1KBB 142,42,43 
43 CU( 2)=RCU(2)*RR*AREM 

CUO ( 2 ) =RCU ( 2 ) 

42 SAUilDCUm 

SAU( 2 1 = C U i 2) 

AY1 i )=ETY { 1)+DE*{AX-ETX(1)) 

TB=A TAN ( DE ) 

XA = A Y ( 1 ) /DX 
JX=XA 
JP= JE+ JX 
AD = J X 
XA=XA— AD 

UP=( l.C-XA)*U( JP)+XA*Ui JP + 1) 

V=SQRT { 2 .C*XQ*ABS ( VA-UP ) ) 

VXt 1 )=V*COSi TB ) 

V Y ( 1)=V*SIN(TB ) 

K = K + 1 
J = 5 1 

DO 118 1 = 1,50 
ETXi J ) = ETX( J— 1 J 
£TY( J ) = E TY ( J— 1 ) 

118 J=J-1 

26 CONTINUE 

21 I F ( IS) 16,16,15 

C ATRAJ CALCULATES THE TRAJECTORIES AFTER THEY HAVE BEEN INIALIZED 

15 CALL TRAJ1K, JE, JED) 

16 K = K + 1 

lF(NTJ-K) 117,116,116 
117 RETURN 

END 
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C AC £LR CALCULATES RHS FOR AX1SYM CASES 

SUBRUU T IKE ACALR (KE.KEOi 

COMMON AREM, AT X { 5 1 ) , A TYI 5 1 ) , AX, AYI51), BASE, CUI51), CUO(51), 

1 DC, DCC , DELY, DX, EPS, ER{40), ET X ( 51 ) , ET Y ( 5 1 ) , H , I AS , J AS ( 20 ) , JD 

2 , JOT , J T ( 3COO ) , KAB, KABB , KAN, KAT, KATT, K AT 1 , KAT2 , KB ( 13) , 

3 KBA, KBB , KBF, KCHI51), KRL, KT(1850), LAST, LB(2), LC<3), MG, 

4 NAJ, NP1T, NRO , NRL , NT, NT J , NURL, NXEP, PTX(99), PTY(99), 

5 RCU<51), RH ( 300C ) , RHDOWN , RHUP , RIN, RX, SAU151), SEM, 

6 SIZE, UI3CCG), UB ( 3000 ) , LRH ( 3000 ) , V A , VAT, VBT , VX(51), 

7 V Y ( 5 1 ) , XD, XEP , XK, XMP, XQM, XR, XT(1850),YEP 
DATA PI/3. 1415527/ 

JE=KE 
J EU=KED 
SWA = 0 .0 

DO 33 J J = 2 , NA J 
J0 = 0 
J=JJ 

C IFCU ( J ) =0 . , NO CONTRIBUTION 
IFICU(J)} 33,33,5 
C IF CUINAJ J , SPECIAL 
5 IFIJ-NAJ) 9,6,33 

6 HT = (RIN-AY{ JJ-1) ) 

XH = AY( J J-l ) 

YL=XH+HT 

JX = AY ( JJ-il/DELY 
NA= JX+JE 
NB=JEO 
KlA = VX ( J J- 1 ) 

WB = WA 
GO TO 23 

9 IF ( KCH ( J ) ) 11,34,34 

C IF TRAJECTORIES HAVE REFLEC TE D , ADD IN MIRROR IMAGE 
34 IFIKCHI J— l)-KCH(J) ) 19,11,15 

11 HT = AY ( J )— AY( J- 1 ) 

IF I HT 1 12,12,13 

12 HT =-HT 
XH=A Y l J ) 

JX=XH/DELY 

NB = AY( J-l l/OELY+l.O 
KA = VX ( J ) 

WB=VXt J-l) 

GO TO 14 

13 XH= A Y l J- 1 3 
JX=Xh/OELY 

R A = V X ( J - 1 ) 

R8 = VX { J ) 

NB=AY( J )/DELY+1.0 

14 NA = JX + JE 
NB = N8 + JE 
YL=XH+HT 
GO TO 23 

IS HT=R IN+RIN-AYl J )-AY ( J-l ) 

J0=1 

R 1 =BASE-AY 1 J ) 

R2 = BASE-AY ( J— 1 ) 
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NB = J ED 

IF { SWA ) 20,20,21 

20 SWA = 1 • C 

FA=R1**2/I Rl**2+R2**2) 

XH= A Y ( J I 
YL = R IN 

JX=Xh/DELY 
NA = J X + J £ 

WA^VXI J) 

WB = VX i J— 1 ) 

GO TO 23 

21 S WA = 0 . 0 

XH = A Y ( J— 1 ) 

YL =R IN 

FA=R2**2/(R1**2+R2**2) 

JX = XH/DE LY 
NA = J X + J E 
WA = VX< J~1 ) 

WB= VX ( J J 

23 DO 32 K=NA, NB 

X A = K— J E 

XUD=(XA-.5)*DELY 
41 XX = XUD +UEL Y 

IF(K-JED) 25,26,26 
26 XX=R IN 

IFiBASE .EQ .RI N ) GO TO 25 
XX= XX + • 5*0 X 
I F { XX. GT. BASE) XX = BASE 
25 YU- AMAX 1 ( XUD , XH ) 

YD = AM IN 1 ( XX , YL ) 

IF ( YD- YU } 32,32,27 
27 XA= .5* ( YD + YU ) 

W=WA4{ XA-XH)*(kB-WA)/HT 
YU= BA S E— YU 
YD = BA SE- YD 
A0 = YU* *2— YD* *2 

AT= 2 . 0* BA S E* ( YL-XH ) +XH ** 2- Y L**2 

XX=BASE-XX 

XUD = B A S £— XUD 

AC=P I*( XUU**2-XX**2) 

IF(JO) 45,45,46 
46 AO=AQ*FA 

45 RH( K i =R H ( K 3 +CU( JJ ) * AO / ( A T*A C*k ) 

32 CONTINUE 
IF(SWA) 33,33,5 

33 CONTINUE 
RETURN 
END 
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APPENDIX D 


FLOW DIAGRAMS 


Flow diagrams (figs. 7 to 28) and brief explanations of all the sections of the program 
and various subroutines, as shown in figure 3, are presented in this appendix. Some of 
the nomenclature used herein is further defined in COMMON STATEMENT SYMBOLS 
(appendix B) and/or in the section Input Data Preparation. 


L1NK0 



Figure 7. - LINKO is used as the Main Program to transfer control between the potential 
field calculation and computation of the ion trajectories and space-charge density. 

) Notation to indicate repetitive operation (FORTRAN DO loop). "(These numbers 
correspond to FORTRAN statement numbers in appendix B.) 


MATRIX 


DO 10 J = 4, M, 3 


DO 11 J = 1, M, 3 



Figure 8. - Subroutine MATRIX solves by Gaussian elimination the matrix equation By » x for the vector y, where B is a tri- 
diagonal matrix of order N stored in an N x 3 array and x is a given N x 1 column vector. 
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CHN14 


DO 1 J = 1, 5 



Figure 9. - Subroutine CHN14 is for data input and initialization. 
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XTCAL 



Figure 10. - Subroutine XTCAL calculates the coefficients in the finite-difference approximation to Poisson's equation. 
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XTCAL 



Figure 10. - Concluded. 
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EVC 


DO 22 KK = 1, 1000 

DO 3 DO 8 (Do a maximum of 1000 iterations) 



Figure 12. - Subroutine EVC calculates the spectral radius of the iteration matrix. The spectral radius is solved 
by a minimax process and as the matrix is two-cyclic, the ratio is calculated every other iteration. 


GUESS 


DO 1 J - 1, NT 



Figure 13. - Subroutine GUESS reads in potential values and assigns 
them to the proper boundaries. 
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CHN12 


Enter 


DO 23 NN = 1, NURL 


1 


Initialize 


Start 
relaxation 
factor at 1 


Set residual 
to zero 
ML = 1 


DO 29 LL = ML, NLIN, 2 

(Iterate by rows, first doing odd, then even) 


Set first 
and last mesh 
point numbers 


DO 38 K = KZ, KC 

(Set up matrix equation for this column) 



Enter 

identity 

U ws 

p. : 

Enter 

coeffic 

i 

matrix 

ients 

elements 




Store 

upper bound 

j 

DO 26 

J = 1, NT 

/Suppress N 
l^space charge^ 


25 




General formula 
to change 
relaxation factor 


Call MATRIX 
to solve for 
potential values 
on this column 


Initialize 

index 



DO 34, 

r K = KZ, KC 



Record total 
suppression 
and print 
bad mesh 
point data 


DO 28 


Yes 


J =2, KNUM 





JN = index 
of test 
point 

\ 

Yes ^ AnyX^ 

■ 

Number of 


test points 

test points 


27 


No 


28 j INIU 



Cali TEST 
for print and 
cycle test 

i 

Continue 

/ 




/ 




/ 


1 


Figure 14. - Subroutine CHN12 calculates the potential field distribution. 
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SUM = length of 
region 

XK = 1 


Eq u ipote ntials\^^J^ 
to be calculated^ * 

|No (SIZE =0) 


^ POTEN No 

out of range " 


Initialize 

potential 

value 

(POTEN) 


heading and ioo*DX 
initialize --r — 


DX - 10 * DX 
SUM = 10* SUM 
RIN = 10 * R IN 
XK = 10*XK 


/RIN <100 \ 
No and ; 
^ \ SUM <100/ 


Yes 

/RIN >100 ’ 
or 

\ SUM >100 


DX = DX/ 10 
SUM = SUM/10 
RIN = RIN/ 10 
XK =XK/10 


DO 22 JJ = 1, JC (Loop on mesh columns) 


1 Yes (L = 1) 

'b\\ values\^_ 
x. printed 1/^” 

|No(L>l) 

DO 24 ,, J = L, 6 

^Set remaining \ 
x-y coordinates 
to zero / 


Print one line of 
x-y coordinates | 


f KS is switch for 
vertical (0) or 
horizontal (1) 
interpolation 

5 l ~7 

AAY is 

y-coordinate on 

mesh column 

KS =0 
JE - JE + 1 

J> No 
(KS - 1) 

x^orizontaHs. 

^ interpolation 


21 j 
Boost JE, 

" JED, AX, JJ 


Yes 

(KS - 0) 


DO 19 K = JE, JED 

f 

J^^rticaT\ 
^ or horizontal 
^\interpolation^ 


J = mesh 
point above 


Horizontal 

s(K = 1) J = ndesh 
— - — ► point to 
^ left 


POTEN 

between U(J) and 
\ U(K) 


Horizontal 
or vertical 
interpolation 


Vertical 
(KS = 0) 

12 


Horizontal 
(KS = 1) 


y-coordinate is AAY 
and interpolate for 
x -coord in ate 


x-coordinate 
* AX + DX and 
interpolate for 
y-coordinate 


DIF = I U(J) -UC(K)| 


U(J) = U(K) 

13 

Store 

coordinate 
of U(K) 


15 <^P rint^> 

16 JYes (L = 6) 
Print out 

one line of — 
x-y coordinates 


Boost 
AAY, K 



L - 

+ 1 




= 0 


Figure 17. - Subroutine EQLINE calculates and prints equipotential lines. 






























CHN13 



Figure 18. - Subroutine CHN13 coordinates the trajectory calculation and space-charge-density-function calculation. 
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TRAJ 



Figure 19. - Subroutine TRAJ calculates the trajectory velocities and coordinates. 


88 



























Mesh 

column in first 
\ electrode/* 


Mesh N. 
[^column in second^ 
^x. electrode/^ 


Calculate index 
of impingement 
coordinate 


Calculate index 
of impingement 
coordinate 


DO 64 K - 1. NTJ 


^Trajectory^x. Yes 
^ initialized/^ * 


/^Trajectory No Store 
terminated or nof/-* trajectory 
\impinging /^ coordinate 


Impingement^^ 
coordinate equal to 
x^width of region/ 


Calculate tube width, 
scale tube current 
and velocities by 
overlap, and reset 
trajectory coordinate 


Set current in 
K plus first 
tube to zero 



Terminate trajectory 
and set tube current 
to zer 


Calculate tube width, scale 


tube current and velocities 


by overlap, and reset 
trajectory coordinate 


Figure 20. - Subroutine CORRCT computes the impingement on the designated electrodes. 








































PEQ 



Figure 21. * Subroutine PEQ calculates the equipotentiai line that is used in conjunction with the emitter to calculate the emitter current density. 
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ARC 



Figure 22. - Subroutine ARC divides the emitter segment into equal increments for intialization of trajectories and calculation of emitter 
current and current density. 


















RUST 



Figure 23. - Subroutine RTEST determines the convergence of the Poisson solution. 
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TROUT 


DO 3 K = 1, NT 



Figure 24. - Subroutine TROUT is used to print RH, the space-charge-density, 
function multiplied by the square of the mesh size divided by 4. 









TRCU 


116 



Figure 25. - Subroutine TRCU initializes the trajectory coordinates and calculates the current in each stream tube for two-dimensional 
problems. 

Subroutine ATRCU initializes the trajectory coordinates and calculates the current in each stream tube for axially symmetric problems. 
The flow chart for subroutine ATRCU is similar to the one for TRCU. The differences are in the internal equations for calculating the 
current, as it is calculated for an annular section in ATRCU but only for a rectangular segment in TRCU, and in the calculation of the 
first stream tube, as it is not possible for it to be symmetric to the upper boundary for axially symmetric problems, while this possibility 
is accounted for in TRCU. 
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TRCU 



Figure 25. - Concluded. 






























































APPENDIX E 


SAMPLE PROBLEMS 

Contained herein is a listing of the necessary input data cards as well as the com- 
puter calculated output data for an axisymmetric problem with current density specified. 
Results are plotted in figure 6 (p. 22). Also included is a two-dimensional problem that 
is solved as a space- charge- limited- flow problem. 
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Input Data Cards 






CAPO COLUMN 

NMVBFRF 




12345673901234567? 

901 

234567 

8 90123456* 7 8 901234^6789 9 12?4 5 678901234667 

390123456780 


7 « 

4 

6 

6 ° 

4 2 

9 1 o 

7 0 

1 

0 A 0 0 1 

6 

1 








CA Q 0 2 

1 

oAM 

PLF DPPBL 

F M Tp 

n r \a n m c j D A 

T F L Fv I o PF 

CFAPCM CFN t 7P 


HFADTNO 


I PM 

T H p ! J C T n P 

PPPPP Av 




hfadt NO 

A X I 

-c YMW 

cio t r 







HFADI no 

CURRENT 

dfnmt 

V cproiFI 

FP 




H f a n t N o 










HFADT N r - 

a 

8 

9 

1 9 






JA 9 










JA.9 

132.91 

1 

p o ^ 


^.26 





AW ,VA 9 

lonn . 

- 

1 Don. 


? 9 9 * 





V/ATVptct 

1.9 

F— 4 1 

.0 r 

-4 

1 .9 

7-4 l.P 

7-4 



PO i 

-1 


5 

-5 

9. 

.25 

.25 

.25 

1 . 

J T = 9 

-1 

1 

5 

~5 

.25 

.25 

.25 

.25 

.75 

J T = 1 4 

“1 

1 

6 

-5 

.25 

.25 

.25 

.75 

. 5 

JT = 1 9 

-1 

2 

c; 

-5 

.25 

.76 

.76 

.76 

.75 

JT = 74 


1 

c 


.25 

.76 

9. 

.26 

9. 

J T = 7 9 

-1 

1 

8 

- 5 

.25 

.26 

.2 6 

.70 

. 5 

J T = o 4 

-1 

1 

5 

-5 

.25 

.25 

.26 

.20 

.25 

J T = 3 9 

-1 

i 

_L 

5 

- 5 

.25 

.25 

9. 

.20 

A 

. 

JT = 44 

-1 

4 

4 

-5 

.25 

.025 

.960 

.25 

.5 

JT = 49 

3 

1 

5 

— 5 

.125 

.25 

.25 

.25 

.25 

J T = 54 

-1 

1 

5 

-1 

.0834 

.25 

.26 

.19 

.5 

JT= 59 

— ] 

j 

5 

-8 

.042 

• 2 ^ 

.76 

* 9 6 

.75 

JT = 64 


1 

F 

— 6 

-.50 

.50 

rs 

• 

9. 

9. 

JT=-69 

-1 

1 

5 

— 5 

.165 

.25 

.25 

.98 

.75 

JT = 74 

.96 

« 

0 5 







ATX 


1 

.0 







AT V 

.625 

• 

40 


.23 





FP 

7 









JA 

1 

2 

0 

4 

-69 

2 1 

3 9 

1 2 


KA* KRF 

1 

1 

n 

1 

74 

2 14 

1 -69 

1 1 1 

64 

0 KA*<R9 

1 

4 

14 

1 

2 





< A = X R 8 

1 

1 

0 

1 

34 

2 19 

1 49 

1 59 6 

19 

0 K A ♦ < n 8 

1 

1 

2 







KA,KP9 

1 

1 

0 

1 

3 9 

2 24 

1 64 

6 74 1 

2 

C K A > K B S 

1 

1 

0 

1 

44 

9 29 

1 7 



K A'> KR S 

1.9 









B A c F 










yp 

1 n o n # 

- 

. 1 O O A # 


9. 





0117 0 0 


co 

CO 
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Output Data Listing 


SAMPLE PROBLEM TO DEMONSTRATE LEWIS RESEARCH CENTER 
ION THRISTCR PROGRAM 
AXI-SYMMETRIC 
CURRENT DENSITY SPECIFIED 


KTiJT) 

KT( JT + U 

KT( JT + 2) 

KT1 JT+3) 

XT ( JT) 

XT ( JT + 1 ) 

XT ( JT+2) 

XT ( JT+3 ) 

XT( JT+4) 

JT 

-1 

1 

5 

“5 

0 . 

0.4827586 

0.2586207 

0.25 86207 

0.2586207 

9 

-1 

1 

5 

-5 

0.2916667 

0.2083333 

0.2500000 

0.2500000 

0.2500000 

14 

-1 

I 

5 

—5 

0.3125000 

0.1875000 

0.2500000 

0.2500000 

0.2500000 

19 

-I 

1 

5 

-5 

0.3750000 

0.1250000 

0.2500000 

0.2500000 

0.2500000 

24 

- I 

1 

5 

-5 

0.6666667 

0 . 

0. 1666667 

0 . 1666667 

0 . 1666667 

29 

-I 

1 

5 

-5 

0.27 77 778 

0.1666667 

0.2469136 

0.3086420 

0.2222222 

34 

-1 

1 

5 

-5 

0.3333333 

0.1111111 

0.24691 36 

0.3086420 

0.2222222 

39 

-1 

1 

5 


0.6153846 

0 . 

0. 1709402 

0.2136752 

0.1538462 

44 

-1 

4 

4 

-5 

0.0651042 

0.2473958 

0.6250000 

0.0625000 

0.C312500 

49 

3 

I 

5 

-5 

0.5797101 

0. 1159420 

0. 1521739 

0. 1521739 

0. 1739130 

54 

”1 

1 

5 

-1 

0.4603994 

0.1063245 

0. 1237932 

0. 3094830 

0.0945295 

59 

-1 

1 

5 

-5 

0.4939706 

0.0672723 

0.0731262 

0. 3656310 

0.0471444 

64 

-4 

1 

5 

-5 

0.5000000 

0.5000000 

0 . 

0 . 

0 . 

-69 

-1 

1 

5 

-5 

0.2225901 

0. 1102924 

0.1617254 

0.5053920 

0. 1098512 

74 


XR= C. 79946046 50 ITERATIONS REQUIRED TC CONVERGE ON XR 

LAPLACE SOLUTION 

AFTER 12 ITERATIONS ON U THE MAXIMUM CHANGE IN U IS 0.00851 VOLTS AND OCCURS AT MESH POINT 12 
60 U VALUES 


WEST POINT NUMBERS POTENTIAL VALUES 


1 

2 

3 

4 

5 

6 

7 

8 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

895. 9544 

766.8770 

9 

10 

1 1 

12 

13 

14 

15 

16 

714.3282 

700.0223 

705.8136 

515.6727 

365.53C3 

292.8042 

273.5507 

2 77. 1696 

17 

18 

19 

20 

21 

22 

23 

24 

38.6665 

-168.9056 

-228.1868 

-229.9403 

-252.4775 

-543.6250 

-918.3443 

-837.2217 

25 

26 

27 

28 

29 

30 

31 

32 

-740.4435 

-771.8125 

- 1000.0000 

-978.6510 

-899.8868 

-863.8362 

- 1000.0000 

-934.5808 

33 

34 

35 

36 

37 

38 

39 

40 

-913.5292 

-862.4311 

-843.0268 

- 1000.0000 

-868.6906 

-797. 9154 

-75 8 . 029 4 

-744.6000 

4 1 

42 

43 

44 

45 

46 

47 

48 

-692.3663 

-658.5857 

-623.7481 

-600.5139 

-592.4560 

- 447 . 7889 

-433.1004 

-423.4583 

49 

50 

51 

52 

53 

54 

55 

56 

-412.1742 

-408.0782 

-221.2972 

-2 18.5140 

-2L3.3258 

-208.9571 

-207.3130 

0 . 

57 

58 

59 

60 

0 

0 

0 

0 

- 0 . 

0 . 

-0 . 

0 . 

0 . 

0 . 

0 . 

0 . 


EQUI POTENT I AL PRINTOUT 


POTENT I AL 

IX, Y) 

1000.0 

I G. 250 , C. ) 

( 0 . , 

0.250) 

( 0 . , 

0.500) 

{ 0 . , 

0.750) 

t 0 . , 

1.00 0 ) 

( 0.250, 

0 . ) 

POTENTIAL 

IX, Yl 

1CC0.Q 

i 0.250, C. ) 

( 0 . 

0 . ) 

{ 0 . , 

0 . > 

10 . , 

0 . ) 

i 0 . 

0 . ) 

t 0 . 

0 . ) 

POTENTIAL 

( X » Y ) 

8 CC.C 

( 0.214, C.500) 

{ 0.175, 

0.750) 

( 0.167, 

1 . 000 ) 

( 0.250, 

0.436 ) 

( 0.420, 

0 . ) 

( 0.313, 

0.2 50 } 

POTENT I AL 

IX, Y) 

600.0 

{ 0.445, 0.250) 

( 0.354, 

0.500) 

( 0.318, 

0.750) 

{ 0.309, 

1 . 000 ) 

( 0.500, 

0.139) 

( 0.562, 

0 . ) 

POTENTIAL 

{ X , Y ) 

400. C 

{ 0.479, C.500) 

I 0.436, 

0.750) 

( 0.426, 

1 . 000 ) 

( 0.500, 

0.44 3 ) 

{ 0.678, 

0 . ) 

( 0.561, 

0.250 ) 

POTENTIAL 

(X,YJ 

200.0 

{ 0.665, C .250) 

{ 0.577, 

0.500) 

( 0.545, 

0. 750 ) 

{ 0.537, 

1 . 000 ) 

( 0.750, 

0.081 ) 

( 0.786, 

0 , ) 

POTENT IAL 

( X , Y ) 

C. 

{ 0.671, C.500) 

{ 0.641, 

0.750 ) 

{ 0.636, 

1 . 000 ) 

( 0.750, 

0.297 ) 

{ 0.881, 

0 . ) 

( 0.767, 

0.250) 

POTENT I AL 

( X , Y ) 

0 . 

i 2.750, C. ) 

I 2.750, 

0.250) 

( 2.750, 

0.500 ) 

( 2. 750, 

0 . 750) 

\ 2.750, 

1 . 000 ) 

i 2.750, 

0 . ) 

POTENT IAL 

t X , Y ) 

c. 

{ 2.750, C.250) 

i 2.750, 

0.500 ) 

( 2.750, 

0.750) 

( 0 . , 

0 . ) 

( 0 . , 

0 . ) 

( 0 . , 

0 . ) 

POTENT IAL 

IX, Y) 

- 200.0 

( 0.736, C . 750) 

i 0.735, 

i.OOO) 

[ 0.750, 

0.631) 

\ 0.975, 

0 . ) 

t 0.852, 

0.25 0) 

t 0.760, 

0 .500) 
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POTENT I AL 

( X , Y ) 

-200.0 

t 2. 524, 

0. ) 

l 2.521, 

0.250) 

( 2.516, 

0.500) 

( 2.511, 

0.750) 

( 2.509, 1.000) 

< 0. 

0. ) 

POTENT I AL 

( X » Y ) 

-400.0 

( 0. 938, 

C .25 0) 

( 0.827, 

0.500) 

{ 0.821, 

0.750) 

{ 0.833, 

1.000) 

( 1.000, 0.127) 

( 1.071, 

0. ) 

POTENTIAL 

(X.Y) 

-4CC.C 

I 2.3C3, 

C. ) 

( 2.293, 

0.250) 

( 2.278, 

0.500) 

( 2.265, 

0.750) 

( 2.260, 1.000) 

( 0. 

0. ) 

POTENT I AL 

IX, Y) 

-600.0 

( 0.894, 

0.500) 

( 0.903, 

0.750) 

( 0.931, 

1.000) 

( 1.000, 

0.288) 

( 1.167, 0. ) 

( 1.031, 

0.250) 

POTENTIAL 

(X.Y) 

-600.0 

( 1.988, 

1.000) 

( 2.000, 

0.766) 

( 2.094, 

0. ) 

( 2.066, 

0.250) 

( 2.030, 0.500) 

( 2.001, 

0.750) 

POTENTIAL 

(X.Y) 

-800.0 

( 0.961, 

0.500) 

( 0.985, 

0.750) 

( 1.000, 

0.421) 

( 1.000, 

0.846) 

( 1.140, 0.250) 

( 1.121, 

1 .000) 

POTENTIAL 

(X.Y ) 

-80C.C 

( 1.250, 

C.031) 

l 1.281, 

0. ) 

( 1.745, 

0.500) 

( 1.649, 

0.750) 

( 1.609, 1.000) 

( 1. 750, 

0.493) 

POTENTIAL 

(X.Y) 

-800.0 

( 1.913, 

0. ) 

( 1.832, 

0.250) 

( 0. 

0. ) 

( 0. , 

0. ) 

( 0. ,0. ) 

( 0. , 

0. ) 

POTENTIAL 

(X.Y) 

-100U.0 

( 1.25C, 

C.250) 

{ 1.250, 

0.250) 

( 1.250, 

0.250) 

( 1.500, 

0. ) 

( 1.250, 0.250) 

( 1.500, 

0. ) 

POTENTIAL 

(X.Y) 

-1000.0 

< 1.750, 

C. ) 

( 1.750, 

0. ) 

( 1.750, 

0. ) 

( 0. , 

0. ) 

( 0. ,0. ) 

< 0. , 

0. ) 


CYCLE 1 


CURRENT DENSITIES ARE CALCULATED USING EOUI POTENT I AL OF 700.02231 VOLTS WHICH HAS X-Y COORDINATES 
i (C. 379.0. 250 2 (0.292,0. 500) 3 (0.258,0.750) 4 (0.250,1.000) 

EMITTER ARC LENGTH , DELTA ARC LENGTH 0.50000 0.12500 

X-Y EMITTER COORDINATES 

1 (0.C50, 0.500 2 (C.C5C, 1.000) 

X-Y EEGIN TRAJ. COORDINATES 


1 ( C. 050, 0.500 2 (0.050,0.625) 3 (0.050,0.750) 4 (0.050,0.875) 


•COORD 

TRAJ NUM 

REFLECTION COUNTER 

Y-COORD 

X-VEL COMP 

Y-VEL COMP 

C.2500 

1 

0 

0.5000 

0.1 8400 E 

05 

0. 



2 

0 

0.6250 

0.19408E 

05 

0 . 



3 

0 

0.7500 

0.20366E 

05 

0 . 



4 

0 

0.8750 

0.20620E 

05 

0 . 


0.5000 

1 

0 

0.5151 

0.30366E 

05 

0 . 294 1 3 E 

04 


2 

0 

0.6337 

0. 312 18 E 

05 

0. 17690 E 

04 


3 

0 

0.7552 

0 . 32045 E 

05 

0. 10844E 

04 


4 

0 

0.8772 

0.32262E 

05 

0.45616E 

03 

C • 7500 

1 

0 

0.5432 

0.41180E 

05 

0.51183E 

04 


2 

0 

0.6501 

0.41707E 

05 

0 • 300 70 E 

04 


3 

0 

0.7648 

0.42218E 

05 

0. 17881 E 

04 


4 

0 

C • 88 09 

0.42239E 

05 

0.65682E 

03 

1.0000 

1 

0 

0.6250 

0. 521 09 E 

05 

0.42150E 

04 


2 

0 

0.6655 

0.51951E 

05 

0.27732E 

04 


3 

0 

0.7731 

0.51534E 

05 

0. 12980E 

04 


4 

0 

0.8825 

0.50922E 

05 

-0. 387 73 E 

02 

1.2500 

1 

0 

0.6425 

0.52960E 

05 

0.31308E 

04 


2 

0 

0.6762 

0 • 52 809 E 

05 

0.17083E 

04 


3 

0 

0.7769 

0.52494E 

05 

0 . 298 59 E 

03 
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4 

0 

0.8802 

0.52253E 

05 

-0.921676 

03 

i.5000 

1 

0 

0.6553 

0. 522 89 E 

05 

0.22822E 

04 


2 

0 

0.6825 

0 . 52 1 83 E 

05 

0.92128E 

03 


3 

0 

0.7769 

0.5I998E 

05 

-0.28952E 

03 


4 

0 

0.8749 

0 . 5 1 8 49 6 

05 

-0.12980E 

04 

1.7500 

1 

0 

0.6650 

0.507616 

05 

0. 16981E 

04 


2 

0 

0.6856 

0.50667E 

05 

0.3 7581 E 

03 


3 

0 

0.7746 

0.50528E 

05 

-0. 685 31 E 

03 


4 

0 

0.8680 

0.50408E 

05 

-0. 15365E 

04 

2.0000 

l 

0 

0.6725 

0.48375E 

05 

0. 12880E 

04 


2 

0 

0.6866 

0.48294E 

05 

-0. 12746E 

02 


3 

0 

0.7704 

0.482I9E 

05 

-0.96813E 

03 


4 

0 

0.8597 

0.48123E 

05 

-0. 170 7 L E 

04 

2.2500 

l 

0 

0.6788 

0.45404E 

05 

0. 105 5 1 E 

04 


2 

0 

0.6859 

0.45326E 

05 

-0.23786E 

03 


3 

0 

0.7647 

0.45299E 

05 

-0. 11364E 

04 


4 

0 

0.8503 

0.45219E 

05 

-0. 18L44E 

04 

2.5000 

I 

0 

0.6845 

0. 4I996E 

05 

0. 9433 IE 

03 


2 

0 

0.6842 

0.41913E 

05 

-0.34820E 

03 


3 

0 

0.7580 

0.41918E 

05 

-0. 1 22 12 E 

04 


4 

0 

0.8397 

0.41846E 

05 

-0. 18713E 

04 

2.7500 

1 

0 

0.6903 

0.38I93E 

05 

0.90961E 

03 


2 

0 

0.6819 

0. 38100E 

05 

-0 • 38200 E 

03 


3 

0 

0.7503 

0.38129E 

05 

-0. 12476E 

04 


4 

0 

0.8280 

0.38059E 

05 

-0. 18897E 

04 


THRUST DISTRIBUTION BY STREAM TUBES (NEWTONS) 


i C. 2 0 • 792058E-06 3 0.128856E-05 4 0.772716E-06 5 0.257337E-06 

TOTAL THRUST (NEWTONS) 0.311067E-05 
TOTAL POWER (WATTS) 0.592785E-01 
INITIAL CURRENT DENSITIES (AMPS/(UNIT H)**2> 
i C. 2 1 • CCOC C06-04 3 L.OOGOOOE-04 4 l. 0000006-04 5 1 . 000000 E- 04 

INITIAL CURRENTS (AMPS) 

10. 2 0. 343612E-04 3 0.2454376-04 4 0.147262E-04 5 0.490874E-05 

TOTAL INITIAL CURRENT =0. 785398E-04 AMPS 

TRANSMITTED CURRENT AT ACCEL. GRID- C.592525E-04 AMPS WHICH IS 75.44 PERCENT UF THE INITIAL CURRENT. 
START OF POISSCN SOLUTION 


HHDUWN* 

0. 

RH ( 9 ) = 6. 7754C11 RHLP- 0. 

U( 

9 ) = 714.3281937 




6 C RH 

VALUES 







I 

MESH 
2 3 

POINT NUMBERS 

4 5 6 7 8 0. 

0. 

RH VALUES 

0 * 0 * 0 . 

0. 

0 . 

4.0842 
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9 

10 

11 

12 

13 

14 15 16 

8.7754 

8.5585 

0 . 0 . 

2.3490 

5. 8381 

5.6656 

0 . 

17 

18 

19 

20 

21 

22 23 24 

0 . 

1.4718 

4.9237 4.6164 

0 . 

0 . 

0 . 

4.6622 

25 

26 

27 

28 

29 

30 31 32 

3.9701 

0 . 

0 . 0 . 

4.5936 

3.7429 

0 . 

0 . 

33 

34 

35 

36 

37 

38 39 40 

0 . 

4.6895 

3.3977 0. 

0 . 

0 . 

4.8723 

3.1376 

41 

42 

43 

44 

45 

46 47 48 

0 . 

0 . 

0. 5.1553 

2.9123 

0 . 

0 . 

0 . 

49 

50 

51 

52 

53 

54 55 56 

5 .5372 

2.7211 

0 . 0 . 

0 . 

6.0326 

2.5649 

0 . 

57 

58 

59 

60 

0 

0 0 0 

0 . 

0 . 

6.6805 2.4476 

0 . 

0 . 

0 . 

0 . 

CYCLE 

2 












RHDGWN* 

0 . 


RH { 

9) 

= 9.C973C64 RHUP 

= 0 . 

U( 9)* 

734.9583206 





CYCLE 

3 












RHDOWN= 

0 . 


RH { 

9 I 

= 9.1C62878 RHUP 

= 0 . 

U I 9 ) = 

735.5056915 





CYCLE 

4 












RFDUWN= 

C. 


RH< 

9) 

= 9.1C65317 RHUP 

= 0 . 

U( 9 ) = 

735.5205154 





CYCLE 

5 












RFDUWN* 

C. 


RHl 

9) 

= 9. 10653 £2 RHUP 

= 0 . 

U ( 9 ) = 

735.5209579 





CYCLE 

6 












RH= AV ER AGE 












RhDOWN = 

0 . 


RHl 

9) 

= 9. 1065348 RHUP' 

= 0 . 

U I 9 ) = 

735.5205307 






CYCLE 7 

THRUST DISTRIBUTION BY STREAM TUBES (NEWTONS) 

1 0. 2 0.6 6 7 437E- 06 3 0.128652E-05 4 0.771840E-06 5 0.257158E-06 

TOTAL THRUST INEWTONS) 0.300296E-05 
TOTAL POWER IWATTS) 0.5715C8E-G1 
INITIAL CURRENT DENSITIES (AMPS/IUNIT H)**2) 

10. 7 1 • OOOOCOE- 04 3 I.CCOOOOE-04 4 1.000000E-04 5 1.000000E-04 

INITIAL CURRENTS {AMPS) 

I 0. 2 0.343612E-C4 3 0.245437E-04 4 0.L47262E-04 5 0.490874E-05 

TOTAL INITIAL CURRENT =0 . 7 8 5 3 S8E- 04 AMPS 

TRANSMITTED CURRENT AT ACCEL. GRID* 0.572759E-04 AMPS WHICH IS 72.93 PERCENT OF THE INITTAL CURRENT. 
RHDUWN* 0. RHl 9>= 9.1065351 RHUP= 0. U( 9)= 735.5207520 

RFOUTPUT IS AVERAGE OF THIS AND PREVIOUS CYCLE 
CONVERGEC POISSCN SOLUTION 
6C RH VALUES 
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MESF 

POINT 

NUMBERS 






RH VALUES 




I 

2 

3 

4 

5 

6 

7 

8 

0 . 

0. 

0 . 

0. 

0 . 

0 . 

0. 

4 . 2 0 7 n 

9 

10 

11 

12 

13 

14 

1 5 

16 

9.1065 

8.9058 

0. 

0. 

2.4018 

5.9097 

5.7551 

0 . 

17 

18 

19 

20 

21 

22 

23 

24 

0 . 

1.5300 

4.8753 

4.5924 

0. 

0. 

0 . 

4.5164 

25 

26 

27 

28 

29 

3C 

31 

32 

3.8570 

0. 

0. 

0. 

4.4666 

3.5396 

0 . 

0. 

33 

34 

35 

36 

37 

38 

39 

40 

0. 

4.56 97 

3.1736 

0. 

0. 

0. 

4. 7608 

2.8615 

41 

42 

4 3 

44 

45 

46 

47 

48 

0. 

0. 

0. 

5,0488 

2.5902 

0 . 

0 . 

0 . 

49 

5C 

51 

52 

53 

54 

55 

56 

5.4304 

2.3564 

0. 

0. 

0. 

5.9130 

2.1576 

0. 

57 

58 

55 

60 

0 

0 

0 

0 

0 . 

0 . 

6.5134 

1.9906 

0 . 

0 . 

0 . 

0 . 

U-F IE L C 

) IS 

AVERAGE OF 

THI S 

AND 

PREVIOUS 

CYCLE 








AFTER 

1 

ITERATIONS 

ON U 

THE 

MAXIMUM CHANGE IN U 

IS 0.00007 VOLTS AND 

OCCURS AT 

MESH POINT 

10 



60 U 

VALUES 
















MESF 

PCINT 

NUMBERS 






POTENTIAL 

VALUFS 




1 

2 

3 

4 

5 

6 

7 

8 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

898. 487 6 

778.7139 

9 

10 

11 

12 

13 

14 

15 

16 

735.5207 

723.6093 

710.8577 

523.2277 

38 1.3708 

319.1874 

303.1791 

285.2541 

17 

18 

19 

20 

21 

22 

23 

24 

47.2805 

-156.2632 

-206.0557 

-204.3045 

-244.7012 

-536.6829 

-917.1026 

-826.3643 

25 

26 

27 

28 

29 

30 

31 

32 

-723 . L 529 

-768.3414 

-1000.0000 

-976.1431 

-884.4876 

-844. 5840 

-1000.0000 

-983.8279 

33 

34 

35 

36 - 

37 

38 

39 

40 

-906.4961 

-842.9459 

-820.5594 - 

1000.0000 

- 864.8654 

-787. 8471 

-735.5510 

-719.6834 

41 

42 

43 

44 

45 

46 

47 

48 

-688.1255 

-652.4277 

-612.1475 

-576.6889 

-566.7832 

-442.8870 

-431. 9096 

-412.6915 

49 

5C 

51 

52 

53 

54 

55 

56 

-389.7871 

-384.6211 

-218.1396 

-214.5997 

-206.3915 

-192.8332 

-191.2206 

0 . 

57 

58 

59 

60 

0 

0 

0 

0 

-0. 

0 . 

- 0 . 

0 . 

0 . 

0 . 

0 . 

0 . 


ECU I POTENT I AL PRINTOUT 


POTENTIAL 

1 X,Y) 

1000.0 

( 0.250, 

C. ) 

( 0. , 

0.250) 

( 0. , 

0.500) 

t 0. , 0.750) 

( 0. , 1.000) 

{ 0.250, 

0. ) 

POTENTIAL 

( X * Y ) 

1C00.C 

I 0.250, 

C. ) 

( 0. 

0. i 

( 0. , 

0. ) 

( 0. ,0. ) 

o 

o 

{ 0. , 

0. ) 

POTENT I AL 

( X , Y ) 

800.0 

i 0.226, 

C.500) 

( 0.189, 0.750) 

( 0.181, 

1.000) 

( 0.250, 0.456) 

( 0.423, 0. ) 

( 0.316, 

0.250) 

POTENT IAL 

f X t Y ) 

600.0 

I 0.449, 

0.250) 

t 0.362, 0.500) 

( 0.331, 

0.750) 

( 0.324, 1.000) 

( 0.500, 0.148) 

{ 0.565, 

0. ) 

POTENTIAL 

IX,Y) 

400.0 

{ 0.488, 

C.500) 

{ 0.451, 0.750) 

( 0.442, 

1.000 ) 

( 0.500, 0.467) 

( 0.683, 0. ) 

( 0.565, 

0.250) 

POTENT IAL 

I X , Y ) 

200.0 

I 0.670, 

C .250) 

I 0.584, 0.500) 

( 0.557, 

0.750) 

I 0.551, 1.000) 

( 0.750, 0.090) 

( 0.790, 

0. ) 

POTENT IAL 

I X,Y> 

C. 

( 0. 677, 

C.500) 

i 0.652, 0.750) 

{ 0.649, 

1.000) 

( 0.750, 0.308) 

{ 0.885, 0. ) 

{ 0.770, 

0.250) 

POTENTIAL 

(X,Y) 

0. 

1 2.75C, 

C. ) 

t 2.750, 0.250) 

( 2.750, 

0.500) 

{ 2.750, 0.750) 

{ 2.750, 1.000) 

( 2.750, 

0. ) 

potential 

( X t Y ) 

0. 

{ 2.750, 

C .250) 

( 2.750, 0.500) 

( 2.750, 

0.750) 

o 

o 

o 

o 

( 0. , 

0. ) 

POTENT IAL 

I X * Y 1 

-200.0 

( 0.747, 

0. 750) 

I 0.748, 1.000) 

( 0.750, 

0.720) 

( 0.979, 0. ) 

[ 0.856, C.250) 

( 0. 764, 

0.500) 

POTENTIAL 

X 

-< 

-200. C 

l 2.491, 

C .750) 

t 2.489, 1.000) 

l 2.500, 

0.618) 

t 2.521, 0. ) 

l 2.517, 0.250) 

( 2.508, 

0.500) 

POTENT IAL 

( X » Y ) 

-400.0 

< 0.941, 

0.250) 

t 0.830, 

0.500) 

( 0.828, 

0.750) 

( 0.844, 1.000) 

{ 1.000, 0.133) 

( 1.074, 

0. ) 

POTENTIAL 

I X * Y ) 

-40G.C 

( 2.236, 

C.750) 

( 2.229, 

o 

o 

o 

• 

t—i 

I 2.250, 

0.639) 

{ 2.298, 0. I 

t 2.287, 0.250) 

( 2.265, 

0.500} 

POTENTIAL 

IX, YJ 

-600.0 

I 0.896, 

C.500) 

I 0.909, 

0.750) 

{ 0.941, 

1.000) 

( 1.000, 0.292) 

o 

o 

I s - 

. 

< 1.034, 

0.250) 

POTENTIAL 

IX, Y) 

-600.0 

I 1.963, 

C.750) 

I 1.946, 1.000) 

i 2.000, 

0.586) 

< 2.090, 0. ) 

( 2.059, 0.250) 

I 2.015, 

0.500) 

PCTENT IAL 

(X,Yi 

-8GC.0 

I 0.962, 

0.500) 

I 0.989, 

0.750) 

( 1.000, 

0.423) 

si- 

CO 

O 

o 

o 

o 

( 1.142, 0.250) 

{ 1.158, 

1.000) 





POTENTIAL 

(X, 

i Y ) 

-8CC.0 

( 1. 

► 250 * 

C.034) 

( 1, 

. 284 , 

0. ) 

( l . 72 A, 

0.500) 

t 1.600, 

0. 750) 

( 1. 

551 , 

1.000) 

( 1. 

750, 

C .461 ) 

POTENTIAL 

(X, 

,Y) 

-800.0 

( 1. 

,910, 

0. > 

( 1. 

.826, 

0.250) 

( 0. 

0. ) 

( 0. 

0. ) 

( 0. 

, 

0. ) 

{ 0. 

, 

0 . ) 

POTENTIAL 

(X, 

► Y) 

-1000.0 

( 1. 

• 25C, 

C.25C) 

( 1, 

.250, 

0.250 ) 

( 1.250, 

0.250) 

{ 1.500, 

0. ) 

( 1. 

250, 

C.250) 

{ 1. 

500, 

0. ) 

POTENTIAL 

(X, 

t Y ) 

-1000.0 

( 1. 

, 750, 

C. ) 

( 1. 

.750, 

0. ) 

( 1.750, 

0. ) 

< 0. , 

0. ) 

( 0. 

, 

0 . ) 

( 0 . 

, 

0 . ) 


CYCLE e 

CURRENT DENSITIES ARE CALCULATED USING E GUI POTENT I AL OF 723.60932 VOLTS rfHlCH HAS X-Y COORDINATES 
1 { 0.489* 0. ) 2 (0.367,0.250) 3 ( 0. 285 , 0 . 500 ) 4 10.257,0.750) 5 10.250,1.000) 

EMITTER ARC LENGTH , DELTA ARC LENGTH 0.50000 0.12500 

X-Y EMITTER COORDINATES 

1 ( 0. C50 , C • 500 1 2 I0.C5C, 1.000) 


X-Y BEGIN TRAJ. COORDINATES 


1 (C.C50 

,0.500) 

2 (0.050,0. 

625) 3 

( 0. C50 , C. 

750) 4 

(0.050,0.875) 


X-CDORD 

TRAJ NUM 

REFLECTION 

COUNTER 

Y-CDURD 

X-VEL COMP 

Y-VEL COMP 

C.2500 

1 

0 


0.5000 

0.17925E 

05 

0. 



2 

0 


0.6250 

C.18779E 

05 

0. 



3 

0 


0.7500 

0. 19596E 

05 

0. 



4 

0 


0.8750 

0.19816E 

05 

0. 


0.5000 

1 

0 


0.5142 

0 . 29983 E 

05 

0.27124E 

04 


2 

C 


0.6327 

0 • 30 721 1 

05 

0. 15246E 

04 


3 

0 


0.7546 

0 • 31441 E 

05 

0 .9 35986 

03 


4 

0 


0.8769 

0.31625E 

05 

0. 39060 E 

03 

C.7500 

1 

0 


0.5405 

0.40959E 

05 

0 . 475 18 E 

04 


2 

0 


0.6470 

0.41405E 

05 

0.25918E 

04 


3 

0 


0.7630 

0.41837E 

05 

0. 15264E 

04 


4 

0 


C. 8800 

0.41827E 

05 

0.52980E 

03 

l.COOO 

1 

0 


0.6250 

0.52064E 

05 

0 • 3 537 1 E 

04 


2 

0 


0.6599 

0. 5 1902E 

05 

0.22244E 

04 


3 

0 


C. 7696 

0.51394E 

05 

0.92984E 

03 


4 

0 


0.8808 

0.50728E 

05 

-0. 25556 E 

03 

1.2500 

1 

0 


0.6389 

0.52874E 

05 

0. 22951 E 

04 


2 

0 


0.6676 

0.52713E 

05 

0. 100 78 E 

04 


3 

0 


0.7714 

0 • 5228 IE 

05 

-0. 1 8571 E 

03 


4 

0 


0.8772 

0. 5 20 1 1 E 

05 

-0. 12283E 

04 

1.5000 

1 

0 


0.6474 

0. 521 33 E 

05 

0. 12675E 

04 


2 

0 


0.6701 

0 . 520 11 E 

05 

0.44516E 

02 


3 

0 


0.7688 

0.51720E 

05 

-0.8923 7 1 

03 


4 

0 


0.8702 

0.51550E 

05 

-0. 16657E 

04 

1.75C0 

1 

0 


0.6517 

0 • 50551 E 

05 

0.50196E 

03 


2 

0 


0.6686 

0.50437E 

05 

-0.67879E 

03 


3 

0 


0.7631 

0.50196E 

05 

-0.14056E 

04 


4 

0 


0.8613 

0. 50059 E 

05 

-0. 19590 E 

04 


o 

05 


2.0000 

1 

0 

0.6527 

0.48U6E 

05 -0.942 7 6 E 

02 


2 

0 

0.6637 

0 • 4801 IE 

05 -0.12500E 

04 


3 

0 

0.7549 

0.47847E 

05 -0.18080E 

04 


4 

0 

C. 8 507 

0 • 47738 E 

05 -0 . 21836 £ 

04 

2.2500 

1 

0 

0.65X1 

0.45129E 

05 -0 . 5 12 20 E 

03 


2 

0 

0.6559 

0 . 450 23 E 

05 -0.16575E 

04 


3 

0 

0.7444 

0.44920E 

05 -0.20944E 

04 


4 

0 

0.8384 

0.44834E 

05 -0.2 341 7 E 

04 

2.50CC 

1 

c 

0.6473 

0.41775E 

05 -0.79243E 

03 


2 

0 

0.6455 

0.41662E 

05 -0.19363E 

04 


3 

0 

0.7318 

0.41608E 

05 -0 . 228 65 E 

04 


4 

0 

0.8246 

0 . 41 543 E 

05 -0.24407E 

04 

2.7500 

l 

0 

0.6420 

0.38169E 

05 -0.90586E 

03 


2 

0 

0.6330 

0.38040E 

05 -0.20509E 

04 


3 

0 

0.7172 

0.38068E 

05 — 0 • 2363 5 E 

04 


4 

0 

C. 8092 

0.3 8033 E 

05 -0.24774E 

04 


THRUST DISTRIBUTION BY STREAM TUBES (NEWTONS) 

10. 2 C.687437E-06 3 0.128652E-05 4 0.771840E-06 5 0.257158E-06 

TOTAL THRUST (NEWTONS) 0.300296E-05 
TOTAL POWER (WATTS) 0.5715C8E-01 
INITIAL CURRENT DENSITIES (AMPS/IUNIT H)**2) 

10. 2 1.0000COE-04 3 1.000000E-04 4 1.000000E-04 5 1.000000E-04 

INITIAL CURRENTS (AMPS) 

1 0. 2 0.343612E-C4 3 0.245437E-04 4 0.147262E-04 5 0.490874E-05 

TOTAL INITIAL CURRENT =0.7853988-04 AMPS 

TRANSMITTED CURRENT AT ACCEL. GR ID= C.572759E-04 AMPS WHICH IS 72.93 PERCENT OF THE INITIAL CURRENT. 


*01* UN I TO 5 » EOF 


FILOOO 


ZECOOO 



Input Data Cards 

card column numbers 

1274567P90l2o4567A9oi?4476" ? Qqni?74767c<9oi?^46E7aQ0i774GA,7£!qn;i?'34667Q90i?74667aQ0 


o 67 4 6 60 4 2 9 1 o 7 o 

5 o 

1 SAMPLF PROBLEM TO DFN n NS T RATF LEWIS Rfaeapch 

CENTER ION THRII7TOR PROGRAM 
TWO-DI MFN A I ONAL 
SPACE CHARGE L I M j T^n 


7 

A 

9 

10 








17?. 91 


inoo. 


.25 







loon. 


-1000, 


200. 







-1 

1 

5 

-5 

0. 


.25 


.25 

.25 

9 

-1 

1 

5 

-5 

.25 


.25 


.25 

.25 

14 

-1 

1 

5 

-5 

.25 


.25 


0. 

.25 

19 

-1 

1 

■ 5 

_C 

.25 


.25 


.25 

.20 

24 

-1 

1 

6 

-5 

.25 


.25 


0. 

.20 

29 

-1 

4 

4 

-6 

.25 


.025 


,060 

.25 

34 

3 

1 

5 

-5 

.125 


.25 


.25 

.25 

39 

-1 

1 

5 

-1 

.0874 


.25 


.27 

.10 

44 

-1 

1 

5 

-5 

.042 


.25 


.25 

.05 

49 

-4 

1 

0 

0 

-.50 


.50 


0. 

0. 

0. -54 

-1 

1 

5 

-5 

.165 


.25 


.25 

.08 

59 

,06 


.05 









.50 


1.0 









.625 


.40 


.23 







9 











1 

5 

0 









1 

1 

0 

1 

59 

2 

24 

1 

29 



2 

1 

-54 

7 

14 

1 

19 





1 

2 

-54 

1 

34 

1 

39 

1 

19 



1 

1 

-54 

1 

1 

1 

44 

1 

14 

1 19 


1 

1 

1 

1 

49 

2 

14 

1 

19 



1 

1 

1 

3 

14 

1 

19 





3 

1 

9 

3 

14 

1 

19 





1 

5 

2 










loon, -lOno. o. 


Output Data Listing 

SAMPLE PROBLEM TO DEMONSTRATE LEV<IS RESEARCH 
CENTER ION THRUSTOR PROGRAM 
TUO-D IMENSIONAL 
SPACE CHARGE LIMITED 


K T ( JT) 

KT( JT+1J 

KTUT + 2) 

KTIJT+3) 

XTt JT) 

XT ( J T + 1 } 

XTI JT+2 > 

XT ( J T+ 3 ) 

XT ( J T + 4 ) 

JT 

-1 

1 

5 

-5 

0. 

0.5000000 

0.2500000 

0. 2500 000 

0.2500000 

6 

-1 

1 

5 

-5 

0.2500000 

0.2500000 

0.2500000 

0. 2500000 

0.2500000 

14 

-1 

1 

5 

-5 

0.5000000 

0. 

0.2500000 

0.2500000 

0.2500000 

1 9 

-1 

1 

5 

-5 

0.2222222 

0.2222222 

0.2469136 

•0.3086420 

0.2222222 

24 

-1 

1 

5 

-5 

0 . 4444444 

0. 

0.2469136 

0.3086420 

0.2222222 

29 

-1 

4 

4 


0.0555556 

0.2777778 

0 • 6060606 

0. 0606061 

0.0333333 

34 

2 

1 

5 


0 . 4444444 

0.2222222 

C. 1666667 

0. 1666 667 

0. 1666667 

39 

-1 

1 

5 

-1 

0.4088604 

0 . 1363958 

0. 1299268 

0.3248L69 

0.0909487 

A4 

- 1 

1 

5 

-5 

0.4653067 

0.0781715 

0.0760870 

0. 3804348 

0.0456522 

49 

-4 

1 

0 

0 

0.5000000 

0.5000000 

0 . 

0 . 

0 . 

-54 

-1 

1 

5 

-6 

0.1967052 

0. 1298254 

0.1632653 

0.5102041 

0. 1077551 

59 


XR = 0.84812816 60 ITERATIONS REQUIRED TO CONVERGE ON XR 

LAPLACE SOLUTION 

AFTER 13 ITERATIONS ON U THE MAXIMUM CHANGE IN U IS 0.00658 VCLTS AND OCCURS AT MESH POINT 45 


6C U VALUES 




MESH 

POINT 

NUMBERS 






POTENTIAL 

VALUES 




1 

2 

3 

4 

5 

6 

7 

8 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

874. 8293 

727.8215 

9 

10 

n 

12 

13 

14 

15 

16 

679.3458 

668.8804 

662. 2943 

449.7594 

2 98.9202 

244.3187 

236.1440 

206.3770 

17 

18 

19 

20 

21 

22 

23 

24 

-37.0053 

-226.2191 

-237.1353 

-212.9401 

-307.4717 

-577.9377 

-929.6562 

-753.7015 

25 

26 

27 

28 

29 

30 

31 

32 

-613.6343 

-788.9688 

-1000.0000 

-953.7076 

-800.2259 

-73 4. 1 93 9 

-1000.0000 

-9 75 ,464 n 

33 

34 

35 

36 

37 

38 

39 

40 

-853.4273 

-759.3019 

-722.6900 

-1000.0000 

-828.1147 

-725.2353 

-660.8605 

-637. 9573 

41 

42 

43 

44 

45 

46 

47 

48 

-659.8769 

-611.7608 

-558.5390 

-520.9495 

-507.4205 

-415.9846 

-400.5105 

-376.2075 

49 

50 

5 i 

52 

53 

54 

55 

56 

-356.9758 

-349.8219 

-203.0415 

-198.0900 

-188. 8058 

-180.9254 

-177.9181 

0. 

57 

58 

59 

60 

0 

0 

C 

0 

-o. 

0. 

-0. 

0. 

0. 

0. 

0. 

0. 


EQUI POTE NT I AL PRINTOUT 


POTENTIAL 

IX, Y) 

1000.0 

I 0.250, 

C. ) 

( 0. , 

0.250) 

( 0. 

0.500) 

( 0. , 

0.750) 

( 0. 

1.000) 1 

1 0.250, 

0. ) 

POTENTIAL 

( X , Y ) 

1000.0 

( 0.250, 

C. ) 

I 0. 

0. ) 

( 0. 

0. ) 

l 0. , 

0. ) 

( 0. , 

0. ) i 

( 0. , 

0. ) 

POTENTIAL 

IX, Y) 

800. C 

( 0.184, 

C. 500) 

( 0.156, 

0.750) 

( 0.151, 

i.OOO) 

{ 0.250, 

0.377) 

I 0.39B, 

0. ) 

I 0.294, 

0.250) 

POTENTIAL 

(X,Y> 

6C0.0 

( 0.412, 

C .250) 

I 0.325, 

0.500) 

( 0.296, 

0.750) 

( 0.290, 

1 .000) 

I 0.500, 

0.073) 

l 0.534, 

0. ) 

POTENTIAL 

I X , Y 1 

4CC.0 

( 0.441, 

C.500) 

( 0.411, 

0.750) 

( 0.405, 

1.000) 

( 0.500, 

0.332) 

I 0.644, 

0. ) 

( 0.526, 

0.250) 

POTENTIAL 

I X , Y ) 

200.0 

( 0. 628, 

C.250) 

( 0.547, 

0.500) 

{ 0.523, 

0.750) 

I 0.520, 

1.000) 

( 0. 750, 

0.007) 

( 0.753, 

0. ) 

POTENTIAL 

( X , Y ) 

c. 

I 0.731, 

C • 250) 

( 0.642, 

0.500) 

( 0.627, 

0.750) 

( 0.631, 

1.000) 

( 0.750, 

0.212) 

I 0. 850, 

0. ) 

POTENT I AL 

( X » Y ) 

0. 

{ 2.750, 

C. ) 

I 2.750, 

0.250) 

l 2.750, 

0.500) 

t 2.750, 

0. 750) 

( 2.750, 

1.000) ' 

( 2. 750, 

n. ) 

POTENTIAL 

I X , Y ) 

0. 

( 2.750, 

C.250) 

< 2.750, 

0.500) 

( 2.750, 

0.750) 

( 0. 

0. ) 

( 0. , 

0. ) 

( 0. , 

0. ) 

POTENTIAL 

( X , Y ) 

-2CC.0 

{ 0. 738, 

C.500) 

( 0.731 , 

0.750) 

I 0.743, 

i.OOO) 

( 0.750, 

0.465) 

I 0.948, 

0. ) 

( 0.825, 

0.250) 

POTENTIAL 

( X , Y ) 

-200.0 

( 2.498, 

C.250) 

( 2.485, 

0.500) 

( 2.473, 

0.750) 

( 2.468, 

1 .000) 

< 2.500, 

0.154) 

( 2.504, 

0. ) 



109 


POTENTIAL 

(X,YJ 

-400.0 

l 0.918, 

C.250) 

( 0.812, 

0.500) 

( 0.829, 

0.750) 

( 0.867, 

1 .000) 

( l.OOO, 

0.086) 

( 1.048, 

0 . ) 

POTENTIAL 

(X,Y) 

-400.0 

( 2.217, 

C.500) 

( 2.184, 

0.750) 

( 2.170, 

1 .000) 

{ 2.250, 

0.255 ) 

( 2. 269, 

0 . ) 

( 2.251, 

0.250) 

POTENTIAL 

( X , Y ) 

-60C.0 

( 0.883, 

C .500 ) 

( 0.926, 

0.750) 

( 0.991, 

1 .000) 

( 1.000, 

0.266) 

( 1.152, 

0 . ) 

( 1.013, 

0.250) 

POTENTIAL 

( X » Y ) 

-600.0 

( 1.938, 

C.500) 

( 1.859, 

0.750) 

( 1.823, 

1 .000) 

( 2.000, 

0.305) 

( 2.061, 

0 . ) 

( 2.014, 

0.250) 

POTENT IAL 

( X , Y ) 

-800.0 

( C.954, 

0.500) 

o 

o 

o 

0.408) 

( 1.000, 

0.684) 

( 1.132, 

0.250) 

( 1.249, 

0.750) 

( 1.250, 

0.013) 

POTENTIAL 

( X , Y ) 

-800.0 

( 1.25C, 

C.751) 

( 1.263 , 

0 . ) 

( 1.251, 

0.750) 

( 1.500, 

0.642) 

( 1.604, 

0.500) 

( 1.750, 

0.318) 

POTENTIAL 

( X,Y) 

-800.0 

( 1.85 7, 

C. ) 

( 1.782, 

0.250 ) 

( 0. 

0 . ) 

( 0. 

0 . ) 

( 0. , 

0 . ) 

( 0 . , 

0 . } 

POTENT IAL 

( X , Y I 

-1000.0 

I 1.250, 

C.250) 

I 1.250, 

0.250) 

1 1.250, 

0.250) 

( 1.500, 

0 . ) 

( 1.250, 

0.250) 

( 1.500, 

0 . ) 

POTENTIAL 

( X , Y ) 

-1000.0 

( 1.75C, 

C. ) 

( 1.750, 

0 . ) 

( 1.750, 

0 . ) 

( 0. 

0 . ) 

( 0. , 

0 . ) 

( 0 . , 

0 , ) 


CYCLE 1 


CURRENT DENSITIES ARE CALCULATED USING E CUI POTENTIAL GF 668.88045 VOLTS WHICH HAS X-Y COORDINATES 
I 10.495,0. ) 2 10.371,0.250) 3 10.284,0.500) 4 10.256,0.750) 5 (0.250,1.000) 

EMITTER ARC LENGTH , DELTA ARC LENGTH 0.50000 0.12500 

X-Y EMITTER COORDINATES 

1 I C.05C, 0. 5CC ) 2 (0.050,1.000) 


X-Y BEGIN TKAJ. COORDINATES 


1 (C.C50,0.5CC) 

2 (0.050,0. 

.625) 3 

(0.050,0. 

750) 4 

(0.' 

050,0.875) 


X-COGRD TRAJ NUM 

REFLECTION 

COUNTER 

Y-COORD 

X-VEL COMP 

Y-VEL COMP 

C.25G0 1 

0 


0.5000 

0.19879E 

05 

0. 


2 

0 


C .625 0 

0 . 20 746 E 

05 

0 . 


3 

0 


0.7500 

0.21577E 

05 

0. 


4 

0 


0.8750 

0.21753E 

05 

0. 


C .5000 1 

0 


0.5132 

0.31909E 

05 

0 . 27388 E 

04 

2 

C 


0.6315 

0 . 32522 E 

05 

0. 13883E 

04 

3 

0 


0. 7537 

0.33124E 

05 

0.79870E 

03 

4 

0 


0.8761 

0.33214E 

05 

0.24523E 

03 

0.75C0 1 

0 


0.5374 

0.42122E 

05 

0 • 441 60 E 

04 

2 

0 


0.6428 

0.42254E 

05 

0. 19884E 

04 

3 

Q 


0.7597 

0.42367E 

05 

0. 10L71E 

04 

4 

0 


0.8772 

0 • 421 71 E 

05 

0 .9442 1 E 

02 

1.C000 1 

0 


0.6250 

0.51713E 

05 

0.18312E 

04 

2 

0 


0.6500 

0.51457E 

05 

0.71830E 

03 

3 

0 


0.7617 

0.50377E 

05 

-0 . 26921 E 

03 

4 

C 


0.8742 

0.49435E 

05 

-0. 12048E 

04 

1.2500 1 

0 


0.6263 

0.52208E 

05 

-0 • 46216 E 

03 

2 

C 


0.6481 

0 . 5 1 983 E 

05 

-0. 151196 

04 

3 

0 


0.7557 

0.51080E 

05 

-0.21590E 

04 

4 

0 


0.8644 

0. 50623 E 

05 

-0.27280E 

04 



o 


1.5000 


1 0 0.6219 0.51217E 05 -0.21946E 04 

2 0 0.6367 0.51046E 05 -0.31806E 04 

3 0 0.7420 0.50486E 05 -0.34097E 04 

4 0 0.8489 0.50214E 05 -0.35370E 04 

1.75C0 1 C 0.6081 0.49568E 05 -0.33575E 04 

2 0 0.6181 0.49412E 05 -0.43143E 04 

3 0 0.7227 0*490 19E 05 -0.42514E 04 

4 0 0.8296 0.48843E 05 -0.40724E 04 

2.C000 1 C 0.5887 0.47225E 05 -0.41589E 04 

2 0 0.5937 0 • 47075 E 05 -0.51043E 04 

3 0 0.6990 0.46836E 05 -0.48469E 04 

4 0 0.8073 0.46757E 05 -0.44556E 04 

2.2500 1 0 0.5647 G.44426E 05 -0.46366E 04 

2 0 0.5644 0.44269E 05 -0.55812E 04 

3 0 0.6713 0.44143E 05 -0.52257E 04 

4 C 0.7821 0 • 441 56 E 05 -0.47086E 04 

2.5000 1 0 0.5369 0.41291E 05 -0.48808E 04 

2 C 0.5310 0.41119E 05 -0.58266E 04 

3 0 0.6401 0.41070E 05 -0.54351E 04 

4 0 0.7541 0 . 41 166fc 05 -0.48533E 04 

2.7500 1 0 0.5059 0.37839E 05 -0.49576E 04 

2 0 0.4938 0.37645E 05 -0.59041E 04 

3 0 0.6053 0 . 37660 E 05 -0.55050E 04 

4 0 0.7233 0. 37833 E 05 -0.49026E 04 

THRUST DISTRIBUTION BY STREAM TUBES (NEWTONS/UNIT H> 

1 C. 2 0.1C1495E-05 3 0.407873E-05 4 0.443956E-05 5 0 . 45 8354E-05 

TUTAL THRUST ( NEWTONS /UN I T H) 0.141168E-04 
TOTAL POWER (WATT/UNIT H) 0.266433E-CC 

INITIAL CURRENT DENSITIES <AMPS/(UNIT H)**2I 

10. 20. 553128E-03 3 0.629138E-03 4 0.6Q3093E-03 5 0.703633E-03 

INITIAL CURRENTS {AMPS/UNIT H) 

1 0. 2 0.651 41 CE- 04 3 0. 786423E-04 4 0.853866E-04 5 0.879541E-04 

TOTAL INITIAL CURRENTS . 32 1 1 24E-03 AMPS/CUNIT H) 

AVERAGE INITIAL CURRENT DENSI TY = 0 .642248E-03 AMPS/IUMT HI **? 

TRANSMITTED CURRENT AT ACCEL. GRID=0. 271 506E-03 AMPS/(UNIT H) WHICH IS 84.55 PERCENT OF THF INITIAL CURRENT. 

START OF POISSON SOLUTION 

RHDQWN= 0. RH { 9)= 54.C54215C RHUP^= 54.0542150 U( 9)= 679.3458099 

6 C RH VALUES 

MESH POINT NUMBERS RH VALUES 

1 2 3 4 5 6 7 8 0. 0. 0. 0. 0. 0. 0. 24.0276 

5 10 11 12 13 14 15 16 54.0542 57.0842 0. 0. 14.3230 35.7465 37.7167 0. 



17 

18 

19 

20 

21 

22 

23 

24 

0. 

9.6175 

29.0589 

29.9902 

0. 

0. 

0. 

25.7305 

25 

26 

21 

28 

29 

30 

31 

32 

24.9610 

0. 

0. 

0. 

26.2311 

22.6091 

0. 

0. 

33 

34 

35 

36 

37 

38 

39 

40 

0.5670 

27.1733 

20.4508 

0. 

0. 

3. 5287 

26.2379 

18.6547 

41 

42 

43 

44 

45 

46 

47 

48 

0. 

0. 

6.4331 

25.8554 

17.2311 

0. 

0. 

10.2188 

49 

50 

51 

52 

53 

54 

55 

56 

25.0918 

16.1353 

0. 

0. 

14.9869 

23.9396 

15.3354 

0. 

57 

58 

59 

60 

0 

0 

0 

0 

0. 

21.0651 

22.2965 

14. 8248 

0. 

0. 

0. 

o. 


CYCLE 2 

RHOOWN= 27.5502651 RH ( 9) = 27. 55C2E51 RHUP= 54.0542150 Ul 9) = 853.3342285 

CYCLE 3 

RF0QWN= 27.5SC2651 RH< 5) = 41.7501826 RHUP = 41.7501826 U( 9) = 758.0848389 

CYCLE 4 

RbC0WN= 3 3 • 77 1CC9C KH ( 9}= 33.771CC9C RHlP = 41.7501826 Ul 9)= 808.6294403 

CYCLE 5 

RHD0WN= 33 . 7 7 ICC 9 C RH 1 9}= 38.C634651 RHLP= 38.0834851 Ul 9)= 779.9215012 

CYCLE 6 
8H = A V fcK A G E 

KHU0WN= 35.5936687 RHI 9)= 36.8386769 RHUP= 38.0834851 Ul 9)= 795.5128098 

CYCLE 7 

THRUST DISTRIBUTION BY STREAM TUBES ( NE WTGNS/UM T H) 

10. 20. 3 0.207486E-05 4 0.225072E-05 5 0.228166E-05 

TOTAL THRUST I NEW TONS/UN I T H) 0.66C724E-05 
TOTAL POWER 1WATT/UNIT H) C.1203C0E-CC 

INITIAL CURRENT DENSITIES (AMPS/IUNIT H)**2) 

10. 2 0.348945E-03 3 0.355756E-03 4 0.358511E-03 5 0.357064E-03 

INITIAL CURRENTS (AMPS/UNIT H) 

1 0. 2 0.436181E— 04 3 0.444695E-04 4 0.448138E-04 5 0.446330E-04 

TOTAL INITIAL CURRENTS. 177534E-03 AMPS/IUNIT H) 

AVERAGE INITIAL CURRENT DENS I T Y = 0. 35 5C69E-03 AMPS/(UMT H)**2 

TRANSMITTED CURRENT AT ACCEL. GRl D=0. 1 31 763E-Q3 AMPS/<UNIT H) WHICH IS 74.22 PERCFNT OF THE INITIAL CURRENT. 
RHDOWN= 35.5938687 RHI 9)= 36.2508059 RHUP= 38.0834851 Ul 9)= 791.0378571 

RHOUTPUT IS AVERAGE OF THIS AND PREVIOUS CYCLE 
CONVERGED PQISSCN SOLUTION 
6C RH VALUES 


1 




MESH 

POINT 1 

NUMBERS 






RH VALUES 




1 

2 

3 

4 

5 

6 

7 

8 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 , 

1 7.9619 

9 

10 

11 

12 

13 

14 

15 

16 

36.5447 

36.5534 

0 . 

0 . 

10.3752 

21.5913 

21.4151 

0 . 

17 

18 

15 

20 

21 

22 

23 

24 

0 . 

7.5518 

16.2042 

15.6423 

0 . 

0 . 

0 . 

13.2213 

25 

26 

27 

28 

29 

30 

31 

32 

12.0545 

0 . 

0 . 

0.9608 

12. 8432 

10.4053 

0 . 

0 . 

33 

34 

35 

36 

37 

38 

35 

40 

2.7482 

12.1621 

8.9666 

0 . 

0 . 

5.3169 

10.9179 

7.8054 

41 

42 

43 

44 

45 

46 

47 

48 

0 . 

0 . 

8.4435 

9.3437 

6.8942 

0 . 

0 . 

11.8661 

49 

50 

51 

52 

53 

54 

55 

56 

7.6627 

6.1796 

0 . 

3.4263 

12.0433 

5.9281 

5.611 3 

0 . 

57 

58 

59 

60 

0 

0 

0 

0 

8.3867 

9.7986 

5.1438 

5.1438 

0 . 

0 . 

0 . 

0 . 

U-FIELD IS 

AVERAGE OF 

THI S 

AND 

PREVIOUS 

CYCLE 








AFTER 

8 

ITERATIONS 

ON U 

THE 

MAXIMUM CHANGE IN U 

IS 0.00473 VOLTS AND 

OCCURS AT 

MESH POINT 

44 



z> 

o 

>0 

VALUES 
















ME S F 

POINT 

NUMBERS 






POTENT IAL 

VALUES 




1 

2 

3 

4 

5 

6 

7 

8 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

1000.0000 

892.071 6 

796.0714 

5 

10 

11 

12 

13 

14 

15 

16 

790.011 7 

790.5113 

696.5851 

501.0987 

395. 551 5 

390.0618 

397.9575 

260.6426 

17 

18 

19 

20 

21 

22 

23 

24 

20.1866 

-146.5260 

-109.6380 

-64.4640 

-257.1412 

-534.4687 

-922.4115 

-682.4409 

25 

26 

27 

28 

29 

30 

31 

32 

-499.1073 

-767.2344 

-1000.0000 

-932.6057 

-705.7480 

-615.2989 

- 1000 . ooon 

-967.7610 

33 

34 

35 

36 

37 

38 

39 

40 

-792.8840 

-644.0203 

-592.2148 - 

1000. OOOC 

-789.0769 

-638.1402 

-533.8799 

-501.3821 

41 

42 

43 

44 

45 

46 

47 

48 

-616.5369 

-550.4070 

-457.9888 

-395.6505 

-376.7771 

-365.3327 

-338.0239 

-281.5294 

45 

50 

51 

52 

53 

54 

55 

56 

-251.3279 

-241.9987 

-168.7470 

-154.8273 

- 126. 2424 

-116.7851 

-113.2813 

0 . 

57 

58 

55 

60 

0 

0 

0 

0 

-0. 

0 . 

- 0 . 

0 . 

0 . 

0 . 

0 . 

0 . 


E GU I POTENTIAL PRINTOUT 


PUTENT IAL 

l X , Y ) 

1000.0 

( 0. 250, 

C. ) 

I 0. , 

0.250) 

( 0. , 

0.500) 

( 0. 

0.750) 

( 0. , 

1.000) 

( 0. 250, 

0. ) 

POTENT IAL 

( X , Y ) 

1000. 0 

( 0. 250* 

C. ) 

< 0. 

0. ) 

l 0. 

0. ) 

1 0. , 

0. ) 

I 0. , 

0. ) 

( 0. , 

0. ) 

POTENT IAL 

t X f Y J 

800.0 

t 0.245, 

C.500) 

I 0.238, 

0.750) 

( 0.239, 

1.000) 

( 0.250, 

0.490) 

( 0.415, 

0. ) 

I 0.309, 

0.250) 

POTENTIAL 

l X * Y ) 

600.0 

{ 0.437, 

C.25C) 

I 0.372, 

0.500) 

I 0.369, 

0. 750) 

( 0.371, 

1.000) 

( 0.500, 

0.124) 

( 0.555, 

0. ) 

POTENTIAL 

( X * Y ) 

400.0 

( 0.497, 

C.500) 

( 0.494, 

0.750) 

( 0.499, 

1.000) 

< 0.500, 

0.489) 

( 0.670, 

0. ) 

( 0.553, 

0.250) 

POTENT IAL 

I X * Y ) 

200.0 

i 0.657, 

C .250 ) 

( C.590, 

0.500) 

( 0.595, 

0.750) 

( 0.607, 

1.000) 

I 0.750, 

0.063 ) 

( 0.779, 

0. ) 

POTENTIAL 

{ X * Y ) 

0. 

t 0.662, 

C.500) 

( 0.695, 

0.750) 

I 0.715, 

1.000) 

{ 0.750, 

0.280) 

{ 0.876, 

0. ) 

i 0.759, 

0.250) 

POTENT IAL 

( X * Y ) 

0. 

( 2.750, 

C. ) 

( 2.750, 

0.250) 

t 2.750, 

0.500) 

I 2.750, 

0. 750) 

( 2.750, 

1.000) 

( 2.750, 

0. ) 

POTENT IAL 

I X, Y) 

0. 

I 2.75C, 

C .250 ) 

I 2.750, 

0.500) 

( 2.750, 

0.750) 

( 0. , 

0. ) 

I 0. , 

0. ) 

( 0. , 

0 . ) 

POTENT IAL 

( X, Y ) 

-200.0 

l 0.972, 

C. ) 

t 0.849, 

0.250) 

< 0.76 7, 

0.500) 

t 0.789, 

0.750) 

( 0.R28, 

1 .000) 

( 2.460, 

0. ) 

POTENTIAL 

(X, Y) 

-200.0 

{ 2.438, 

C .25 0) 

< 2.381, 

0.500) 

( 2. 345, 

0.750) 

( 2. 332, 

1.000) 

( 0. 

0. ) 

i 0. , 

0. ) 

POTENT IAL 

I X * Y ) 

-40C.C 

I 0.939, 

C.250) 

I 0.832, 

0.500) 

{ 0.877, 

0. 1 50) 

( 0.943, 

1.000) 

I 1.000, 

0.129) 

{ 1.070, 

0. ) 

POTENT IAL 

{ X * Y ) 

-400.0 

{ 1. 992, 

0. 750) 

( 1.953, 

1.000) 

( 2.000, 

0.733) 

( 2.215, 

0. ) 

( 2.L77, 

0.250) 

( 2. 082 , 

0,500) 

POTENT IAL 

IX, Y) 

-600.0 

( 0. 896, 

C.500) 

l 0.964, 

0.750) 

( 1.000, 

0.292) 

I 1.000, 

0.862 ) 

( 1.168, 

0. ) 

( 1.035 , 

0.250) 

POTENTIAL 

( X,Y) 

-60C.0 

{ 1.217, 

l.OCC) 

( 1.416, 

. 

c 

o 

o 

I 1.500, 

0.962 ) 

( 1.600, 

0. 750) 

( 1.750, 

0.591 ) 

( 1.948, 

0.250) 

POTENTIAL 

IX, Y) 

-600.0 

I 1.803, 

C.500) 

o 

o 

o 

r\j 

0.063) 

I 2.016, 

0. ) 

( 0. , 

0. ) 

I 0. , 

0. ) 

( 0. , 

C. ) 



POTENTIAL 

( X , Y ) 

-800.0 

( C* 961 1 

C.500) 

( 1.000, 

0.421 ) 

( 1.000, 

0.628) 

( 1.143, 

0.250) 

( 1.250, 

0.035) 

{ 1.250, 

0.646) 

POTENTIAL 

( X , Y ) 

-800.0 

1 1. 285, 

0. ) 

( 1.487, 

0.500) 

( 1.500, 

0.490) 

{ 1.735, 

0.250) 

( 1.750, 

0.237) 

( 1.880, 

0. ) 

POTENTIAL 

( X , Y ) 

-1000.0 

( 1.25C, 

C.250) 

( 1.250, 

0.250) 

( 1.250, 

0.250) 

( 1.500, 

0. ) 

( 1.250, 

C.250) 

( 1.500, 

0. ) 

POTENTIAL 

IX, Y) 

1 

r— 

o 

o 

o 

• 

o 

( 1.750, 

0. ) 

( 1.750, 

0. ) 

( 1.750, 

0. ) 

{ 0. 

0. ) 

( 0. , 

0. ) 

( 0. 

0. ) 


CYCLE 8 


CURRENT DENSITIES ARE CALCULATED USING ECUi POTENT IAL OF 790.51134 VOLTS WHICH HAS X-Y COORDINATES 
l (0.423*0. ) 2 (0. 315,0. 25Ci 3 (0.253,0.500) 4 (0.249,0.750) 5 (0.250,1.000) 

EMITTER ARC LENGTH , DELTA ARC LENGTH 0.50000 0.12500 

X-Y EMITTER COORDINATES 

I 1 C. 050, C . 50C 1 2 (0.050,1.000) 


X-Y BEGIN TRAJ. COORDINATES 


1 1 C.Q50 

,0.500 

2 (0.050,0. 

,625) 3 

(0.050,0. 

750) 4 

(O.i 

050,0.875 ) 


X-COGRC 

TRAJ NUM 

REFLECTION 

COUNTER 

Y-COORD 

X-VEL COMP 

Y-VEL COMP 

C • 2500 

1 

0 


C. 5000 

0.17207E 

05 

0 . 



2 

0 


0.6250 

0.17335E 

05 

0 . 



3 

0 


0.7500 

0.1 746 1 E 

05 

0 . 



4 

0 


0.8750 

0. 1745 1 E 

05 

0 . 


0.5000 

1 

0 


0.5085 

0.29625E 

05 

0.15987E 

04 


2 

0 


0.6259 

0.29692E 

05 

0. 177 52 E 

03 


3 

0 


0.7501 

0.29759E 

05 

0.24140E 

02 


4 

0 


0.8743 

0.29663E 

05 

-0.128806 

03 

0.7500 

1 

0 


0.5235 

0.40752E 

05 

0.26086E 

04 


2 

0 


0.626C 

0.40468E 

05 

-0.14826E 

03 


3 

0 


0.7487 

0 . 40 1 41 E' 

05 

-0.41414E 

03 


4 

0 


0.8714 

0.39735E 

05 

-0.68234E 

03 

1.0000 

1 

0 


- 1.0000 

0 . 52451 E 

05 

0. 279846 

04 


2 

0 


0.6250 

0 • 5 1 i 04 E 

05 

-0.22442E 

04 


3 

0 


0.7407 

0* 49490E 

05 

-0.24776E 

04 


4 

0 


0.8621 

0.48134E 

05 

-0.25935E 

04 

1.2500 

1 

0 


-1.0000 

0 • 524 51 E 

05 

0. 27984E 

04 


2 

0 


0.6062 

0.51334E 

05 

-0.54675E 

04 


3 

0 


0.7210 

0.49816E 

05 

-0.531576 

04 


4 

0 


0.8431 

0.49051E 

05 

-0.47907E 

04 

1.5000 

1 

0 


- 1.0000 

0 • 52451 fc 

05 

0.27984E 

04 


2 

0 


0.5726 

0.49728E 

05 

-0.81195E 

04 


3 

0 


0.6885 

0. 4 8698 E 

05 

-0.75167E 

04 


4 

0 


0.8150 

0.48313E 

05 

-0.61550E 

04 

1.7500 

1 

0 


- 1.0000 

0. 52451 E 

05 

0.27984E 

04 


2 

0 


0.5255 

0.47549E 

05 

-0.10188E 

05 


3 

0 


0.6447 

0.46797E 

05 

-0.92106E 

04 


4 

0 


0.7799 

0.46686E 

05 

-0.71860E 

04 


2.C000 

1 

0 

-i.OOGO 

0. 52451E 

05 

0. 27984E 

04 


2 

0 

0 • 4-659 

0.44683E 

05 

-0.U796E 

05 


3 

0 

0* 5905 

0.44238E 

05 

-0.10527E 

05 


4 

0 

0.7382 

0. 444 8 1 1 

05 

-0.80020E 

04 

2 . 2 50C 

1 

0 

-1.0000 

0. 52451 & 

05 

0.27984E 

04 


2 

0 

0.3942 

0.41543 E 

05 

-0. 12941 E 

05 


3 

0 

0.5264 

0.41373E 

05 

-0.114356 

05 


4 

c 

0.6903 

0.41980E 

05 

-0.85679E 

04 

2 . 500C 

1 

0 

-1.0000 

0.52451E 

05 

0.27984E 

04 


2 

0 

0.3110 

0.38425E 

05 

-0. 1 3661 E 

05 


3 

0 

0.4530 

0 . 38523 E 

05 

-0.12009E 

05 


4 

0 

0.6368 

0.39454E 

05 

-0.88704E 

04 

2 ml 500 

1 

0 

-1.0000 

0 * 52451 E 

05 

0.279846 

04 


2 

0 

0.2177 

0.35432E 

05 

-0. 13905E 

05 


Q 

0 

0.3715 

0.35863E 

05 

-0.12230E 

05 


4 

0 

0.5786 

0.37138E 

05 

-0.89573E 

04 


THRUST DISTRIBUTION BY STREAM TUBES < NEWTONS/UM T HI 

10. 20. 3 0.209802E-05 4 0.226902E-05 5 0 . 2 30 1 77E- 09 

TOTAL THRUST ( NEW TONS /UN I T H) 0.666880E-05 
TOTAL POWER (WATT/UNIT H) 0.121546E-0C 

INITIAL CURRENT DENSITIES (AMPS/IUNIT H)**2) 

1 0. 2 C.350769E-03 3 0.357920E-03 4 0.361040E-03 5 0.359966E-03 

INITIAL CURRENTS IAMPS/UNIT H) 

1 0. 2 0.438462E-04 3 0.447401E-04 4 0.451300E-04 3 0.449937E-04 

TOTAL INITIAL CURRENTS . I787I2E-03 AMPS/IUNIT H) 

AVERAGE INITIAL CURRENT DEN SI TY = 0. 3 57424E-03 AMPS/IUMT H)**2 


TRANSMITTED CURRENT AT ACCEL. GR 1 0 = 0 . 1 32 853E-03 AMPS/<UNIT H) WHICH IS 


74.34 PERCENT Of THE INITIAL CURRENT 
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